#include "ContigGraphAnalyzer.h"


ContigGraphAnalyzer::ContigGraphAnalyzer(void)
{
	m_contigs = new map<int, Contig*>();
	m_headList = new map<string, list<int>* >;	
	m_tailList = new map<string, list<int>* >;
	m_headHashMap = new hash_map<const char*, list<int>*, hash<const char*>, eqName>;
	m_tailHashMap = new hash_map<const char*, list<int>*, hash<const char*>, eqName>;
	m_contigsName = new map<string, Contig*>;
	m_fileType = SOAP;
	m_totalGaps = 0;
	m_multiPath = 0;
	m_covThreshold = 0;
	m_outputPath = "results/";
	m_aimLinearValue = new vector<int>;
	m_aimQuadValue = new vector< vector<int>* >;
	m_simplePath = "";
}


ContigGraphAnalyzer::~ContigGraphAnalyzer(void)
{
	map<int, Contig*>::iterator mapIter = m_contigs->begin();
	while( mapIter != m_contigs->end() ){
		delete mapIter->second;
		m_contigs->erase( mapIter++ );
	}
	delete m_contigs;

	m_headHashMap->clear();
	delete m_headHashMap;
	m_tailHashMap->clear();
	delete m_tailHashMap;

	m_headList->clear();
	delete m_headList;
	m_tailList->clear();
	delete m_tailList;

	m_contigsName->clear();
	delete m_contigsName;

	m_aimLinearValue->clear();
	delete m_aimLinearValue;

	for( vector<vector<int>*>::iterator quadIter = m_aimQuadValue->begin(); quadIter != m_aimQuadValue->end(); quadIter++ )
	{
		(*quadIter)->clear();
		delete *quadIter;
	}
	m_aimQuadValue->clear();
	delete m_aimQuadValue;
}

// set coverage
void ContigGraphAnalyzer::SetCovThreshold( double cov )
{
	m_covThreshold = cov;
}

// analyze graph file
int ContigGraphAnalyzer::ReadVelvetGraphFile( string graphFileName )
{
	ifstream graphReader( graphFileName.c_str() );
	m_avgCov = 0;

	if( graphReader == NULL )
		return -1;

	// analyze the contig number and kmer size
	string line;
	getline( graphReader, line );
	vector<string> content;
	Split( line, "\t", &content );

	m_numOfContig = atoi( content.at( 0 ).c_str() );
	m_sumLength = atof( content.at( 1 ).c_str() );
	Configure::OVERLAP = atoi( content.at( 2 ).c_str() ) - 1;

	string plusSeq = "";
	string minusSeq = "";
	while( getline( graphReader, line ) ){
		if( line.substr( 0, 4 ) == "NODE" ){
			// get contig information
			getline( graphReader, plusSeq );
			getline( graphReader, minusSeq );
			GetVelvetContigInfo( line, plusSeq, minusSeq );
			
		}
		else if( line.substr( 0, 3 ) == "ARC" ){
			// get edge information
			GetVelvetEdgeInfo( line );
		}
	}

	m_avgCov /= m_sumLength;

	return 0;
}

// get contig information
void ContigGraphAnalyzer::GetVelvetContigInfo( string line, string plusSeq, string minusSeq )
{	
	// read contig ID and length
	vector<string> content;
	Split( line, "\t", &content );
	Contig *newContig = new Contig( atoi( content.at( 1 ).c_str() ), atoi( content.at( 2 ).c_str() ) + Configure::OVERLAP );
	newContig->SetPlusSeq( plusSeq );
	newContig->SetMinusSeq( minusSeq );
	m_contigs->insert( pair<int, Contig*>( atoi( content.at( 1 ).c_str() ), newContig ) );

	// calculate coverage
	//double cov = ( atof( content.at( 3 ).c_str() ) + atof( content.at( 5 ).c_str() ) ) / newContig->GetLength();
	double cov = atof( content.at( 3 ).c_str() ) / atoi( content.at( 2 ).c_str() );
	newContig->SetCov( cov );
	m_avgCov += cov * atoi( content.at( 2 ).c_str() );
}

// get edge information
void ContigGraphAnalyzer::GetVelvetEdgeInfo( string line )
{
	// read contig ID 
	vector<string> content;
	Split( line, "\t", &content );

	int startContigID = atoi( content.at( 1 ).c_str() );
	int endContigID = atoi( content.at( 2 ).c_str() );
	int size = atoi( content.at( 3 ).c_str() );
	Contig *endContig = (*m_contigs)[ abs( endContigID ) ];
	Contig *startContig = (*m_contigs)[ abs( startContigID ) ];

	if( startContigID > 0 ){
		// add right edge to start contig
		startContig->AddRightEdge( endContigID, endContig, size );
	}
	else{
		// add left edge to start contig
		startContig->AddLeftEdge( -endContigID, endContig, size );
	}

	if( endContigID > 0 ){
		// add left edge to end contig
		endContig->AddLeftEdge( startContigID, startContig, size );
	}
	else{
		// add right edge to end contig
		endContig->AddRightEdge( -startContigID, startContig, size );
	}
}

// find interesting contigs
void ContigGraphAnalyzer::FindRepeats(){
	ofstream contigWriter( "possibleRepeats.txt" );

	string tempString = "Total contigs are: " + itos( m_numOfContig ) + "\n";
	contigWriter.write( tempString.c_str(), tempString.length() );

	tempString = "Average coverage is " + itos( m_avgCov ) + "\n";
	contigWriter.write( tempString.c_str(), tempString.length() );

	if( contigWriter == NULL ){
		cout<<"error writing possibleRepeats.txt"<<endl;
		return;
	}

	tempString = "\n\nContigs with high coverage and multiple neighbors\n";
	contigWriter.write( tempString.c_str(), tempString.length() );
	map<int, Contig*>::iterator mapIter = m_contigs->begin();
	while( mapIter != m_contigs->end() ){
		if( mapIter->second->GetCov() > m_avgCov * 1.5 && mapIter->second->HasMultipleNeighbors() ){
			tempString = "";
			tempString += "ID: " + itos( mapIter->second->GetID() ) + "\t";
			tempString += "Length: " + itos( mapIter->second->GetLength() ) + "\t";
			tempString += "Cov: " + itos( mapIter->second->GetCov() ) + "\n";
			contigWriter.write( tempString.c_str(), tempString.length() );
		}
		mapIter++;
	}

	tempString = "\n\nContigs with low coverage but multiple neighbors\n";
	contigWriter.write( tempString.c_str(), tempString.length() );
	mapIter = m_contigs->begin();
	while( mapIter != m_contigs->end() ){
		if( mapIter->second->GetCov() <= m_avgCov * 1.5 && mapIter->second->HasMultipleNeighbors() ){
			tempString = "";
			tempString += "ID: " + itos( mapIter->second->GetID() ) + "\t";
			tempString += "Length: " + itos( mapIter->second->GetLength() ) + "\t";
			tempString += "Cov: " + itos( mapIter->second->GetCov() ) + "\n";
			contigWriter.write( tempString.c_str(), tempString.length() );
		}
		mapIter++;
	}

	tempString = "\n\nContigs with high coverage but no multiple neighbors\n";
	contigWriter.write( tempString.c_str(), tempString.length() );
	mapIter = m_contigs->begin();
	while( mapIter != m_contigs->end() ){
		if( mapIter->second->GetCov() > m_avgCov * 1.5 && !mapIter->second->HasMultipleNeighbors() ){
			tempString = "";
			tempString += "ID: " + itos( mapIter->second->GetID() ) + "\t";
			tempString += "Length: " + itos( mapIter->second->GetLength() ) + "\t";
			tempString += "Cov: " + itos( mapIter->second->GetCov() ) + "\t";
			tempString += "number of neighbors: " + itos( mapIter->second->GetNumberOfNeighbors() ) + "\n";
			contigWriter.write( tempString.c_str(), tempString.length() );
		}
		mapIter++;
	}

	contigWriter.close();
}

// analyze SOAPdenovo graph file
int ContigGraphAnalyzer::ReadSoapFile( string graphFileName, string contigFileName )
{
	// read contig file to get contig ID and length
	ifstream contigReader( contigFileName.c_str() );
	m_avgCov = 0;
	m_sumLength = 0;

	if( contigReader == NULL ){
		cout<<"Cannot open file: "<<contigFileName<<endl;
		return -1;
	}

	string line;
	int numOfContig = 0;
	while( getline( contigReader, line ) ){
		if( line.substr( 0, 1 ) == ">" ){
			char nameChar[ 100 ];
			char lengthChar[ 50 ];
			int length = 0;
			sscanf( line.c_str(), "%s %s %d", nameChar, lengthChar, &length );

			string name = string( nameChar );
			name = name.substr( 1, name.length() - 1 );
			name = "C" + name;
			vector<string> *covVec = new vector<string>;
			Split( line, "_", covVec );
			double cov = atof( covVec->at( 1 ).c_str() );
			covVec->clear();
			delete covVec;

			// read contig ID and length
			Contig *newContig = new Contig( numOfContig, length );
			newContig->SetName( name );
			newContig->SetCov( cov );
			m_contigs->insert( pair<int, Contig*>( numOfContig, newContig ) );
			m_contigsName->insert( pair<string, Contig*>( name, newContig ) );
			/*if( cov > 0 ){
				cout<<line<<endl;
				cout<<name<<"\t"<<length<<"\tcov: "<<cov<<endl;
			}*/
			
			// calculate coverage
			m_avgCov += cov * length;
			numOfContig++;
			m_sumLength += length;
		}
	}

	cout<<"There are "<<numOfContig<<" contigs\n";

	m_avgCov /= m_sumLength;
	contigReader.close();
	//cout<<"average cov: "<<m_avgCov<<endl;

	ReadSoapGraphFile( graphFileName );
}

// analyze SOAPdenovo graph file
int ContigGraphAnalyzer::ReadSoapGraphFile( string graphFileName ){
	// read updated.edge file to get the overlap information
	ifstream graphReader( graphFileName.c_str() );

	if( graphReader == NULL ){
		cout<<"Cannot open file: "<<graphFileName<<endl;
		return -1;
	}

	// find the first valid line
	string line;
	char buffer[ 20 ];
	int firstLength = (*m_contigs)[ 0 ]->GetLength() - Configure::OVERLAP;
	sprintf( buffer, "%d,", firstLength );
	string firstLine = string( buffer );

	while( getline( graphReader, line ) ){
		if( line.at( 0 ) == '>' ){
			// check the length
			char first[ 100 ];
			char second[ 100 ];
			sscanf( line.c_str(), "%s %s", first, second );
			string secondString = string( second );
			if( secondString.substr( 0, firstLine.length() ) == firstLine ){
				break;
			}
		}
	}

	//cout<<line<<endl;
	// read each line and save the first and last kmer
	int contigID = 0;
	string firstKmer, lastKmer;
	do{
		GetTwoKmer( line, firstKmer, lastKmer );

		//cout<<line<<endl;
		//cout<<firstKmer<<"\t"<<lastKmer<<endl;

		// get length
		int spacePos = line.find( " " );
		int commaPos = line.find( "," );
		string lengthString = line.substr( spacePos + 1, commaPos - spacePos - 1 );
		int length = atoi( lengthString.c_str() );
		//cout<<"length : "<<length<<endl;
		if( length + Configure::OVERLAP != (*m_contigs)[ contigID ]->GetLength() ){
			cout<<"Find unmatch!"<<endl;
			cout<<"current contigID is: "<<contigID<<endl;
			cout<<(length+Configure::OVERLAP)<<"\t"<<(*m_contigs)[ contigID ]->GetLength()<<endl;
			cout<<"contig name is "<<(*m_contigs)[ contigID ]->GetName()<<endl;
			cout<<line<<endl;
			break;
		}
		
		// save head and tail for PLUS orientation
		InsertKmerIntoMap( firstKmer, m_headList, contigID, PLUS );
		InsertKmerIntoMap( lastKmer, m_tailList, contigID, PLUS );
		//InsertKmerIntoHashMap( firstKmer, m_headHashMap, contigID, PLUS );
		//InsertKmerIntoHashMap( lastKmer, m_tailHashMap, contigID, PLUS );
		//cout<<firstKmer<<endl;
		/*if( m_headHashMap->find( firstKmer.c_str() ) == m_headHashMap->end() ){
			cout<<"Error: not find kmer: "<<firstKmer<<endl;
		}*/

		// check if it is palindrome
		if( line.find( ",0," ) == string::npos ){
			// not palindrom, deal with next line
			getline( graphReader, line );
			GetTwoKmer( line, firstKmer, lastKmer );
			// save head and tail
			InsertKmerIntoMap( firstKmer, m_headList, contigID, MINUS );
			InsertKmerIntoMap( lastKmer, m_tailList, contigID, MINUS );
			//InsertKmerIntoHashMap( firstKmer, m_headHashMap, contigID, MINUS );
			//InsertKmerIntoHashMap( lastKmer, m_tailHashMap, contigID, MINUS );
		}

		contigID++;
	}while( getline( graphReader, line ) );

	// find all overlap information
	FindOverlapOfSOAP();

	graphReader.close();
	

	return 0;
}

// analyze ABySS file
int ContigGraphAnalyzer::ReadAbyssFile( string graphFileName, string contigFileName ){
	// read contig file to get contig ID and length
	ifstream contigReader( contigFileName.c_str() );
	m_avgCov = 0;
	m_sumLength = 0;

	if( contigReader == NULL ){
		cout<<"Cannot open file: "<<contigFileName<<endl;
		return -1;
	}

	string line;
	int numOfContig = 0;
	while( getline( contigReader, line ) ){
		if( line.substr( 0, 1 ) == ">" ){
			char nameChar[ 100 ];
			int length = 0;
			int kmer = 0;
			sscanf( line.c_str(), "%s %d %d", nameChar, &length, &kmer );
			cout<<"length: "<<length<<" kmer: "<<kmer<<endl;

			string name = string( nameChar );
			name = name.substr( 1, name.length() - 1 );
			// convert total kmer to kmer coverage
			double cov = kmer / (length - Configure::OVERLAP);
			cout<<"cov: "<<cov<<endl;

			// read contig ID and length
			Contig *newContig = new Contig( numOfContig, length );
			newContig->SetName( name );
			newContig->SetCov( cov );
			m_contigs->insert( pair<int, Contig*>( numOfContig, newContig ) );
			m_contigsName->insert( pair<string, Contig*>( name, newContig ) );
			
			// calculate coverage
			m_avgCov += cov * length;
			numOfContig++;
			m_sumLength += length;
		}
	}

	cout<<"There are "<<numOfContig<<" contigs\n";

	m_avgCov /= m_sumLength;
	contigReader.close();

	ReadAbyssGraphFile( graphFileName );

	return 0;
}

// analyze ABySS graph file
int ContigGraphAnalyzer::ReadAbyssGraphFile( string graphFileName ){
	// read dot file to get the overlap information
	ifstream graphReader( graphFileName.c_str() );

	if( graphReader == NULL ){
		cout<<"Cannot open file: "<<graphFileName<<endl;
		return -1;
	}

	// handle each 
	string line;

	while( getline( graphReader, line ) ){
		// split the line
		//string firstContigName, arrow, secondContigName, dis;
		vector<string> *edgeLine = new vector<string>;
		Split( line, " ", edgeLine );
		if( edgeLine->at( 0 ) == "edge" )
		{ // get kmer value
			Configure::OVERLAP = atoi( edgeLine->at( 1 ).substr( 4, edgeLine->at( 1 ).length() - 5 ).c_str() ) + 1;
		}
		else if( edgeLine->size() >= 3 ){
			if( edgeLine->at( 1 ) == "->" )
			{ // valid line
				cout<<line<<endl;

				string firstContigOri = edgeLine->at( 0 ).substr( edgeLine->at( 0 ).length() - 2, 1 );
				Contig *firstContig = (*m_contigsName)[ edgeLine->at( 0 ).substr( 1, edgeLine->at( 0 ).length() - 3 ) ];

				string secondContigOri = edgeLine->at( 2 ).substr( edgeLine->at( 2 ).length() - 2, 1 );
				Contig *secondContig = (*m_contigsName)[ edgeLine->at( 2 ).substr( 1, edgeLine->at( 2 ).length() - 3 ) ];

				int dis = -( Configure::OVERLAP - 1 );
				if( edgeLine->size() > 3 )
				{
					dis = -atoi( edgeLine->at( 3 ).substr( 4, edgeLine->at( 3 ).length() - 5 ).c_str() );
				}

				if( firstContig->GetID() <= secondContig->GetID() )
				{ // create the edge, and avoid duplicating edges
					if( firstContigOri == "+" )
					{ // create right edge for first Contig
						if( secondContigOri == "+" ){
							firstContig->AddRightEdge( secondContig->GetID(), secondContig, dis );
							// create left edge for second contig
							secondContig->AddLeftEdge( firstContig->GetID(), firstContig, dis );
						}
						else{
							firstContig->AddRightEdge( -secondContig->GetID(), secondContig, dis );
							// create right edge for second contig
							secondContig->AddRightEdge( -firstContig->GetID(), firstContig, dis );
						}
					}
					else{
						// create left edge for first contig
						if( secondContigOri == "+" ){
							firstContig->AddLeftEdge( -secondContig->GetID(), secondContig, dis );
							// add left edge for second contig
							secondContig->AddLeftEdge( -firstContig->GetID(), firstContig, dis );
						}
						else{
							firstContig->AddLeftEdge( secondContig->GetID(), secondContig, dis );
							// add right edge for second contig
							secondContig->AddRightEdge( firstContig->GetID(), firstContig, dis );
						}

					}
				}
			}
		}
		edgeLine->clear();
		delete edgeLine;
	}

	graphReader.close();

	return 0;
}

// get two kmers from the line
void ContigGraphAnalyzer::GetTwoKmer( string line, string &firstKmer, string &lastKmer ){
	vector<string> *temp = new vector<string>;
	Split( line, ",", temp );
	lastKmer = (*temp)[ 3 ];

	// get the first kmer, by removing the first space
	int spacePos = (*temp)[ 2 ].find( " " );
	firstKmer = (*temp)[ 2 ].substr( spacePos + 1, (*temp)[ 2 ].length() );
	//cout<<firstKmer<<"\t"<<lastKmer<<endl;

	temp->clear();
	delete temp;
}

// insert the kmer into map
void ContigGraphAnalyzer::InsertKmerIntoMap( string kmer, map<string, list<int>* >* kmerList, int contigID, int ori ){
	// save head and tail
	map<string, list<int>* >::iterator firstIter = kmerList->find( kmer );
	if( firstIter == kmerList->end() ){
		// not exist
		list<int> *newList = new list<int>;
		newList->push_back( contigID );
		newList->push_back( ori );
		kmerList->insert( pair<string, list<int>* >(kmer, newList) );
	}
	else{
		firstIter->second->push_back( contigID );
		firstIter->second->push_back( ori );
	}
}

void ContigGraphAnalyzer::InsertKmerIntoHashMap( string kmer, hash_map<const char*, list<int>*, hash<const char*>, eqName> *kmerHashMap, int contigID, int ori ){
	// save head and tail
	hash_map<const char*, list<int>*, hash<const char*>, eqName>::iterator firstIter = kmerHashMap->find( kmer.c_str() );
	if( firstIter == kmerHashMap->end() ){
		// not exist
		list<int> *newList = new list<int>;
		newList->push_back( contigID );
		newList->push_back( ori );
		kmerHashMap->insert( pair<const char*, list<int>* >(kmer.c_str(), newList) );
		
	}
	else{
		firstIter->second->push_back( contigID );
		firstIter->second->push_back( ori );
	}
}

// find all overlap information of SOAPdenovo
void ContigGraphAnalyzer::FindOverlapOfSOAP(){
	/*int count = 0;
	cout<<"total first kmer: "<<m_headHashMap->size()<<endl;
	cout<<"total last kmer: "<<m_tailHashMap->size()<<endl;
	hash_map<const char*, list<int>*, hash<const char*>, eqName>::iterator headIter = m_headHashMap->begin();
	for( headIter = m_headHashMap->begin(); headIter != m_headHashMap->end(); headIter++ ){
		string kmer = string( (*headIter).first );
		count++;
		hash_map<const char*, list<int>*, hash<const char*>, eqName>::iterator tailIter = m_tailHashMap->find( (*headIter).first );
		if( tailIter != m_tailHashMap->end() ){
			// there is overlap
			cout<<"kmer is: "<<kmer<<endl;
			cout<<"Tail: ";
			list<int>::iterator listIter = (*tailIter).second->begin();
			while( listIter != (*tailIter).second->end() ){
				cout<<(*m_contigs)[ *listIter ]->GetName()<<" ";
				listIter++;
				cout<<(*listIter)<<"\t";
				listIter++;
			}
			cout<<endl;
			cout<<"Head: ";
			listIter = (*headIter).second->begin();
			while( listIter != (*headIter).second->end() ){
				cout<<(*m_contigs)[ *listIter ]->GetName()<<" ";
				listIter++;
				cout<<(*listIter)<<"\t";
				listIter++;
			}
			cout<<endl<<endl;

			if( count > 3 )
				break;
		}
	}*/

	cout<<"total first kmer: "<<m_headList->size()<<endl;
	cout<<"total last kmer: "<<m_tailList->size()<<endl;

	int numOfEdge = 0;

	map<string, list<int>* >::iterator headIter = m_headList->begin();
	map<string, list<int>*>::iterator tailIter;
	for( headIter = m_headList->begin(); headIter != m_headList->end(); headIter++ ){
		string kmer = (*headIter).first;
		tailIter = m_tailList->find( (*headIter).first );
		if( tailIter != m_tailList->end() ){
			// there is overlap
			//cout<<"kmer is: "<<kmer<<endl;
			//cout<<"Tail: ";
			list<int>::iterator headListIter = (*headIter).second->begin();
			while( headListIter != (*headIter).second->end() ){
				Contig *headContig = (*m_contigs)[ *headListIter ];
				headListIter++;
				int headOri = *headListIter;
				list<int>::iterator tailListIter = (*tailIter).second->begin();
				while( tailListIter != (*tailIter).second->end() ){
					Contig *tailContig = (*m_contigs)[ *tailListIter ];
					tailListIter++;
					int tailOri = *tailListIter;
					if( tailContig->GetID() < headContig->GetID() ){
						// add edge
						numOfEdge++;
						if( headOri == PLUS ){
							if( tailOri == PLUS ){
								// + +
								headContig->AddLeftEdge( tailContig->GetID(), tailContig, 1 );
								tailContig->AddRightEdge( headContig->GetID(), headContig, 1 );
							}
							else{
								// - +
								headContig->AddLeftEdge( -tailContig->GetID(), tailContig, 1 );
								tailContig->AddLeftEdge( -headContig->GetID(), headContig, 1 );
							}
						}
						else{
							if( tailOri == PLUS ){
								// + -
								headContig->AddRightEdge( -tailContig->GetID(), tailContig, 1 );
								tailContig->AddRightEdge( -headContig->GetID(), headContig, 1 );
							}
							else{
								// - -
								headContig->AddRightEdge( tailContig->GetID(), tailContig, 1 );
								tailContig->AddLeftEdge( headContig->GetID(), headContig, 1 );
							}
						}
					}
					
					tailListIter++;
				}
				
				headListIter++;
			}
		}
	}

	cout<<"There are "<<numOfEdge<<" edges\n";

	// delete all kmer map
	headIter = m_headList->begin();
	while( headIter != m_headList->end() ){
		(*headIter).second->clear();
		delete (*headIter).second;
		m_headList->erase( headIter++ );
	}

	tailIter = m_tailList->begin();
	while( tailIter != m_tailList->end() ){
		(*tailIter).second->clear();
		delete (*tailIter).second;
		m_tailList->erase( tailIter++ );
	}

	//PrintEdges();
}

// print edges of all contigs
void ContigGraphAnalyzer::PrintEdges()
{
	int numOfContig = m_contigs->size();
	for( int i = 0; i < numOfContig; i++ ){
		cout<<"Contig "<<(*m_contigs)[ i ]->GetName()<<endl;
		Contig *currentContig = (*m_contigs)[ i ];
		cout<<"left edges:"<<endl;
		PrintEdges( currentContig->GetEdge( LEFT ) );
		cout<<"right edges:"<<endl;
		PrintEdges( currentContig->GetEdge( RIGHT ) );
		cout<<endl;
	}
}

// print certain edges
void ContigGraphAnalyzer::PrintEdges( map<int, pair<Contig*, int> > *edges )
{
	map<int, pair<Contig*, int> >::iterator iter = edges->begin();
	while( iter != edges->end() ){
		cout<<"\t"<<(*iter).second.first->GetName()<<"\t";
		if( (*iter).first < 0 )
			cout<<"MINUS\t";
		else
			cout<<"PLUS\t";
		iter++;
	}
	cout<<endl;
}

// print dot file of current graph
void ContigGraphAnalyzer::PrintDotFile()
{
	// find subgraph to print
	set<Contig*> *visitedContigs = new set<Contig*>;
	set<Contig*> *subgraph = new set<Contig*>;
	list<Contig*> *possibleContigs = new list<Contig*>;
	int subgraphID = 1;

	int numOfContig = 0;
	for( map<int, Contig*>::iterator contigIter = m_contigs->begin(); contigIter != m_contigs->end(); contigIter++ ){
		Contig *currentContig = contigIter->second;

		// do not consider the contigs with low coverage
		if( currentContig->GetCov() < m_covThreshold )
			continue;

		// if visited before, do not use this contig
		if( visitedContigs->find( currentContig ) != visitedContigs->end() )
			continue;

		subgraph->clear();
		possibleContigs->clear();
		possibleContigs->push_back( currentContig );
		visitedContigs->insert( currentContig );
		numOfContig = 1;
		while( !possibleContigs->empty() ){
			numOfContig++;
			//cout<<"current number of contig: "<<numOfContig<<endl;
			//cout<<"visited contigs size is "<<visitedContigs->size()<<endl;

			Contig *c = possibleContigs->front();
			possibleContigs->pop_front();

			subgraph->insert( c );
			map<int, pair<Contig*, int> > *edges = c->GetEdge( LEFT );
			for( map<int, pair<Contig*, int> >::iterator edgeIter = edges->begin(); edgeIter != edges->end(); edgeIter++ ){
				Contig *otherContig = edgeIter->second.first;
				// do not consider the contigs with low coverage
				if( otherContig->GetCov() < m_covThreshold )
					continue;

				if( visitedContigs->find( otherContig ) == visitedContigs->end() ){
					// add this contig
					possibleContigs->push_back( otherContig );
					visitedContigs->insert( otherContig );
				}
			}
			edges = c->GetEdge( RIGHT );
			for( map<int, pair<Contig*, int> >::iterator edgeIter = edges->begin(); edgeIter != edges->end(); edgeIter++ ){
				Contig *otherContig = edgeIter->second.first;
				// do not consider the contigs with low coverage
				if( otherContig->GetCov() < m_covThreshold )
					continue;

				if( visitedContigs->find( otherContig ) == visitedContigs->end() ){
					// add this contig
					possibleContigs->push_back( otherContig );
					visitedContigs->insert( otherContig );
				}
			}
		}

		// find a subgraph, output it if it is not a single node
		if( subgraph->size() > 1 )
			PrintSubgraph( subgraphID++, subgraph, "subgraph" );
	}
	
	visitedContigs->clear();
	delete visitedContigs;
	subgraph->clear();
	delete subgraph;
	possibleContigs->clear();
	delete possibleContigs;
}

// print subgraph
int ContigGraphAnalyzer::PrintSubgraph( int id, set<Contig*> *graph, string name )
{
	string fileName = m_outputPath + name + itos( id ) + ".dot";
	FILE *graphFile = fopen( fileName.c_str(), "w" );
	fprintf( graphFile, "graph %s%d {\nrankdir=LR;\noverlap=vpsc\n", name.c_str(), id );

	// print contigs
	for( set<Contig*>::iterator contigIter = graph->begin(); contigIter != graph->end(); contigIter++ ){
		Contig *currentContig = (*contigIter);
		fprintf( graphFile, "\"%s length: %d cov: %.1f\"", currentContig->GetName().c_str(), currentContig->GetLength(), currentContig->GetCov() );
		if( currentContig->IfInScaffold() || currentContig->IfVec() )
			fprintf( graphFile, " [color = red, style=filled];\n" );
		else
			fprintf( graphFile, ";\n" );
	}

	// print edges
	int numOfEdge = 0;
	for( set<Contig*>::iterator contigIter = graph->begin(); contigIter != graph->end(); contigIter++ ){
		Contig *currentContig = (*contigIter);
		numOfEdge += PrintEdges( graphFile, currentContig, currentContig->GetEdge( LEFT ), LEFT, graph );
		numOfEdge += PrintEdges( graphFile, currentContig, currentContig->GetEdge( RIGHT ), RIGHT, graph );
	}

	fprintf( graphFile, "}\n" );
	fclose( graphFile );

	if( numOfEdge == 0 )
		m_disconnectedGap++;

	return numOfEdge;
}

// print edges to file
// only print edges whose other contig ID is bigger than id of c
int ContigGraphAnalyzer::PrintEdges( FILE *file, Contig *c, map<int, pair<Contig*, int> > *edges, int type, set<Contig*> *subgraph )
{
	string arrowtail = "";
	if( type == LEFT )
		arrowtail = "vee";
	else
		arrowtail = "crow";

	int numOfEdge = 0;
	for( map<int, pair<Contig*, int> >::iterator iter = edges->begin(); iter != edges->end(); iter++ ){
		if( c->GetID() > iter->second.first->GetID() )
			continue;

		if( type == RIGHT && c->GetID() == abs( iter->first ) )
			continue;

		Contig *otherContig = iter->second.first;

		// do not output the contigs which is not in this subgraph
		if( subgraph->find( otherContig ) == subgraph->end() )
			continue;

		// do not output contigs with low coverage
		if( otherContig->GetCov() < m_covThreshold )
			continue;

		string arrowhead = "";
		if( type == LEFT ){
			if( iter->first > 0 )
				arrowhead = "crow";
			else
				arrowhead = "vee";
		}
		else{
			if( iter->first > 0 )
				arrowhead = "vee";
			else
				arrowhead = "crow";
		}

		fprintf( file, "\"%s length: %d cov: %.1f\" -- \"%s length: %d cov: %.1f\" ", c->GetName().c_str(), c->GetLength(), c->GetCov(), 
			otherContig->GetName().c_str(), otherContig->GetLength(), otherContig->GetCov() );
		fprintf( file, "[fontsize=30, arrowsize=1, len=3.1, penwidth=1.80, arrowtail=%s, arrowhead=%s, style=dashed, dir=both, color=purple, label = \"%d\"];\n", 
			arrowtail.c_str(), arrowhead.c_str(), iter->second.second );

		/*if( c->IfVec() && otherContig->IfVec() && CheckIfNormal( otherContig, c, arrowhead, arrowtail ) )
		{ // normal edge
			fprintf( file, "[fontsize=30, arrowsize=1, len=3.1, penwidth=2.80, arrowtail=%s, arrowhead=%s, style=solid, dir=both, color=black];\n", 
				arrowtail.c_str(), arrowhead.c_str());
		}
		else
			fprintf( file, "[fontsize=30, arrowsize=1, len=3.1, penwidth=1.80, arrowtail=%s, arrowhead=%s, style=dashed, dir=both, color=purple];\n", 
				arrowtail.c_str(), arrowhead.c_str());*/

		numOfEdge++;
	}

	return numOfEdge;
}

// vector related funcitons
bool ContigGraphAnalyzer::CheckIfNormal( Contig *head, Contig *tail, string arrowHead, string arrowTail )
{
	int firstPos = 0;
	int secondPos = 0;

	if( arrowHead == "vee" ){
		if( arrowTail == "vee" ){
			// <-  ->
			if( head->GetVecOri() != tail->GetVecOri() ){
				if( head->GetVecOri() == PLUS ){
					secondPos = head->GetStartInVec();
					firstPos = tail->GetEndInVec();
				}
				else{
					firstPos = head->GetEndInVec();
					secondPos = tail->GetStartInVec();
				}
			}
		}
		else{
			// <-  <-
			if( head->GetVecOri() == tail->GetVecOri() ){
				if( head->GetVecOri() == PLUS ){
					secondPos = head->GetStartInVec();
					firstPos = tail->GetEndInVec();
				}
				else{
					firstPos = head->GetEndInVec();
					secondPos = tail->GetStartInVec();
				}	
			}
		}
	}
	else{
		if( arrowTail == "vee" ){
			// <- <-
			if( head->GetVecOri() == tail->GetVecOri() ){
				if( head->GetVecOri() == PLUS ){
					firstPos = head->GetEndInVec();
					secondPos = tail->GetStartInVec();
				}
				else{
					secondPos = head->GetStartInVec();
					firstPos = tail->GetEndInVec();
				}					
			}
		}
		else{
			// ->  <-
			if( head->GetVecOri() != tail->GetVecOri() ){
				if( head->GetVecOri() == PLUS ){
					firstPos = head->GetEndInVec();
					secondPos = tail->GetStartInVec();
				}
				else{
					secondPos = head->GetStartInVec();
					firstPos = tail->GetEndInVec();
				}
			}
		}
	}

	if( firstPos - secondPos + 1 == Configure::OVERLAP )
		return true;
	else
		return false;
}

// read scaffold file
void ContigGraphAnalyzer::ReadScaffoldFile( string fileName )
{
	ifstream scafReader( fileName.c_str() );
	m_totalFilledGaps = 0;

	if( scafReader == NULL )
		return;

	ofstream scafWriter( (fileName + ".gapFilled").c_str() );

	// read all scaffold to calculate the occurence for all used contigs
	string line;
	int numOfContig = 0;
	Contig *firstContig;
	while( getline( scafReader, line ) ){
		if( line.substr( 0, 1 ) == ">" ){
			numOfContig == 0;
		}
		else{
			vector<string> *temp = new vector<string>;
			Split( line, "\t", temp );
			Contig *currentContig;
			if( m_fileType == SOAP )
				currentContig = (*m_contigsName)[ (*temp)[ 0 ] ];
			else{
				vector<string> *tempName = new vector<string>;
				Split( (*temp)[ 0 ], "_", tempName );
				currentContig = (*m_contigs)[ atoi( (*tempName)[ 1 ].c_str() ) ];
				tempName->clear();
				delete tempName;
			}
			temp->clear();
			delete temp;

			numOfContig++;
			if( numOfContig == 1 )
				firstContig = currentContig;
			if( numOfContig == 2 )
				firstContig->AddOccu();

			currentContig->AddOccu();
		}
	}
	scafReader.close();

	scafReader.open( fileName.c_str() );

	// analyze the contig number and kmer size
	string preline, nextline;;
	vector<string> *content = new vector<string>;
	vector<string> *scafLine = new vector<string>;
	string gapString;
	int scafID = 0;
	list<Gap*> *oriScaf = new list<Gap*>;		// the original scaffold
	list<Gap*> *sortGap = new list<Gap*>;		// the sorted gaps
	string scafName;
	string outputString;

	while( getline( scafReader, nextline ) ){
		if( nextline.substr( 0, 1 ) == ">" ){
			if( scafID > 0 ){
				// calculate gap for previous scaffold
				sortGap->sort( small_gap );

				// fill gap
				//cout<<"scaffold\n";
				for( list<Gap*>::iterator iter = sortGap->begin(); iter != sortGap->end(); iter++ )
				{
					//cout<<(*iter)->GetGapSize()<<endl;
					FillGap( *iter );
				}
				//cout<<endl;

				// output and delete previouse gap
				scafWriter.write( scafName.c_str(), scafName.length() );
				Split( (*oriScaf->begin())->GetStartString(), "\t", scafLine );
				outputString = (*scafLine)[ 0 ] + "\t" + (*scafLine)[ 1 ] + "\t" + (*scafLine)[ 2 ] + "\t";
				gapString = (*scafLine)[ 3 ];
				scafWriter.write( outputString.c_str(), outputString.length() );
				for( list<Gap*>::iterator iter = oriScaf->begin(); iter != oriScaf->end(); iter++ ){
					if( (*iter)->GetGapString() != "" ){
						outputString = "0\n";
						if( (*iter)->GetGapString() != "Connect" )
							outputString += (*iter)->GetGapString();
					}
					else
						outputString = gapString + "\n";
					Split( (*iter)->GetEndString(), "\t", scafLine );
					outputString += (*scafLine)[ 0 ] + "\t" + (*scafLine)[ 1 ] + "\t" + (*scafLine)[ 2 ] + "\t";
					gapString = (*scafLine)[ 3 ];
					scafWriter.write( outputString.c_str(), outputString.length() );
					delete *iter;
				}
				// output the end of scaffolds
				outputString = "0\n";
				scafWriter.write( outputString.c_str(), outputString.length() );
			}

			// start a new scaffold
			scafID++;
			oriScaf->clear();
			sortGap->clear();
			scafName = nextline + "\n";
			getline( scafReader, preline );
		}
		else{
			// generate a gap
			Gap *newGap = GenerateGap( preline, nextline );

			oriScaf->push_back( newGap );
			sortGap->push_back( newGap );
			preline = nextline;
		}
	}

	// calculate gap for last scaffold
	sortGap->sort( small_gap );
	// fill gap
	//cout<<"scaffold\n";
	for( list<Gap*>::iterator iter = sortGap->begin(); iter != sortGap->end(); iter++ )
	{
		//cout<<(*iter)->GetGapSize()<<endl;
		FillGap( *iter );
	}
	//cout<<endl;

	// delete previouse gap
	scafWriter.write( scafName.c_str(), scafName.length() );
	Split( (*oriScaf->begin())->GetStartString(), "\t", scafLine );
	outputString = (*scafLine)[ 0 ] + "\t" + (*scafLine)[ 1 ] + "\t" + (*scafLine)[ 2 ] + "\t";
	gapString = (*scafLine)[ 3 ];
	scafWriter.write( outputString.c_str(), outputString.length() );
	for( list<Gap*>::iterator iter = oriScaf->begin(); iter != oriScaf->end(); iter++ ){
		if( (*iter)->GetGapString() != "" ){
			outputString = "0\n";
			if( (*iter)->GetGapString() != "Connect" )
				outputString += (*iter)->GetGapString();
		}
		else
			outputString = gapString + "\n";
		Split( (*iter)->GetEndString(), "\t", scafLine );
		outputString += (*scafLine)[ 0 ] + "\t" + (*scafLine)[ 1 ] + "\t" + (*scafLine)[ 2 ] + "\t";
		gapString = (*scafLine)[ 3 ];
		scafWriter.write( outputString.c_str(), outputString.length() );
		delete *iter;
	}

	// output the end of scaffolds
	outputString = "0\n";
	scafWriter.write( outputString.c_str(), outputString.length() );

	cout<<"There are "<<m_totalGaps<<" gaps in total\n";
	cout<<m_totalFilledGaps<<" gaps have been filled\n";
	cout<<m_multiPath<<" gaps have multi paths\n";

	scafLine->clear();
	delete scafLine;
	content->clear();
	delete content;
	oriScaf->clear();
	delete oriScaf;
	sortGap->clear();
	delete sortGap;

	scafWriter.close();
}

// generate Gap
Gap* ContigGraphAnalyzer::GenerateGap( string preline, string nextline )
{
	vector<string> *preTemp = new vector<string>;
	vector<string> *nextTemp = new vector<string>;
	Split( preline, "\t", preTemp );
	Split( nextline, "\t", nextTemp );

	// get first contig
	Contig *firstContig;
	vector<string> *temp = new vector<string>;
	if( m_fileType == VELVET ){
		Split( (*preTemp)[ 0 ], "_", temp );
		firstContig = (*m_contigs)[ atoi( (*temp)[ 1 ].c_str() ) ];
	}
	else{
		firstContig = (*m_contigsName)[ (*preTemp)[ 0 ] ];
	}
	int firstOri = PLUS;
	if( (*preTemp)[ 1 ] == "EB" )
		firstOri = MINUS;

	int gap = atoi( (*preTemp)[ 3 ].c_str() );

	// get second contig
	Contig *secondContig;
	if( m_fileType == VELVET ){
		temp->clear();
		Split( (*nextTemp)[ 0 ], "_", temp );
		secondContig = (*m_contigs)[ atoi( (*temp)[ 1 ].c_str() ) ];
	}
	else{
		secondContig = (*m_contigsName)[ (*nextTemp)[ 0 ] ];
	}
	int secondOri = PLUS;
	if( (*nextTemp)[ 1 ] == "EB" )
		secondOri = MINUS;

	Gap *newGap = new Gap( firstContig, firstOri, secondContig, secondOri, gap, "" );
	newGap->SetEndString( nextline );
	newGap->SetStartString( preline );

	temp->clear();
	delete temp;
	preTemp->clear();
	delete preTemp;
	nextTemp->clear();
	delete nextTemp;

	return newGap;
}

// calculate average coverage
void ContigGraphAnalyzer::CalculateAvgCov()
{
	double avgCov = 0;
	double sumLength = 0;
	for( map<int, Contig*>::iterator iter = m_contigs->begin(); iter != m_contigs->end(); iter++ ){
		//if( (*iter).second->GetLength() > 1000 ){
			avgCov += iter->second->GetCov() * (double) iter->second->GetLength();
			sumLength += iter->second->GetLength();
		//}
	}
	avgCov = avgCov / sumLength;
	//cout<<"average cov is "<<avgCov<<endl;

	for( map<int, Contig*>::iterator iter = m_contigs->begin(); iter != m_contigs->end(); iter++ ){
		int copyNum = (int) ( iter->second->GetCov() / avgCov );
		if( copyNum < 1 )
			copyNum = 1;

		copyNum++;

		iter->second->SetCopyNum( copyNum );
	}
}

// fill in gaps
void ContigGraphAnalyzer::FillGap( Gap *gap )
{
	Contig *startContig = gap->GetStartContig();
	Contig *endContig = gap->GetEndContig();
	cout<<"start contig: "<<startContig->GetName()<<endl;
	cout<<"end contig: "<<endContig->GetName()<<endl;

	list<string> *paths = new list<string>;
	PathNode *newNode = new PathNode( startContig, gap->GetStartOri(), 0, 0, "" );

	cout<<"current gap size: "<<gap->GetGapSize()<<endl;
	FindPath( newNode, paths, gap, true );

	if( paths->empty() || ( Configure::STRINGENT && paths->size() > 1 ) ){
		cout<<"Cannot fill in this gap"<<endl;
	}
	else{
		cout<<"Find "<<paths->size()<<" paths to fill in this gap\n";
		m_totalFilledGaps++;
		SeletePath( gap, paths );
		GenerateGapSeq( gap );
		cout<<gap->GetGapString()<<endl<<endl;
	}

	cout<<"finish current gap\n";
	paths->clear();
	delete paths;
	delete newNode;
}

void ContigGraphAnalyzer::GenerateGapSeq( Gap *g )
{
	string gapScaf = g->GetGapString();
	vector<string> *lines = new vector<string>;
	Split( gapScaf, "\n", lines );
	string scafSeq = "";
	for( int i = 0; i < lines->size(); i++ )
	{
		vector<string> *oneLine = new vector<string>;
		Split( lines->at( i ), "\t", oneLine );
		vector<string> *names = new vector<string>;
		Split( oneLine->at( 0 ), "_", names );
		Contig *c = (*m_contigs)[ atoi( names->at( 1 ).c_str() ) ];
		if( oneLine->at( 1 ) == "BE" )
			scafSeq += c->GetPlusSeq();
		else
			scafSeq += c->GetMinusSeq();
		oneLine->clear();
		delete oneLine;
		names->clear();
		delete names;
	}

	if( scafSeq.length() > Configure::OVERLAP ){
		scafSeq = scafSeq.substr( 0, scafSeq.length() - Configure::OVERLAP );
		scafSeq += "\n";
		g->SetGapSeq( scafSeq );
		g->SetCalculatedGapSize( scafSeq.length() );
	}
	else
	{
		g->SetCalculatedGapSize( scafSeq.length() - Configure::OVERLAP );
		scafSeq = "\n";
		g->SetGapSeq( scafSeq );
	}
	lines->clear();
	delete lines;
}

// find path in iterative way
// if exit due to multiple paths, return false; else return true;
bool ContigGraphAnalyzer::FindPath( PathNode *node, list<string> *paths, Gap *gap, bool isStart )
{
	// if the distance is to long, do not consider this node
	if( node->GetDis() > gap->GetMaxGapSize() ){ 
		if( paths->empty() && isStart )
		{
			m_disconnectedGap++;
		}

		return true;
	}

	map<int, pair<Contig*, int> > *edges;
	if( node->GetOri() == PLUS )
		edges = node->GetContig()->GetEdge( RIGHT );
	else
		edges = node->GetContig()->GetEdge( LEFT );

	for( map<int, pair<Contig*, int> >::iterator edgeIter = edges->begin(); edgeIter != edges->end(); edgeIter++ ){
		Contig *newContig = edgeIter->second.first;
		if( gap->GetSubgraph()->find( newContig ) == gap->GetSubgraph()->end() && newContig != gap->GetEndContig()  )
		{ // this contig is not within this subgraph
			continue;
		}

		int newContigOri = edgeIter->first;
		if( node->GetOri() == MINUS )
			newContigOri = -newContigOri;

		if( newContigOri > 0 )
			newContigOri = PLUS;
		else
			newContigOri = MINUS;

		if( newContig == gap->GetEndContig() && newContigOri == gap->GetEndOri() ){
			// find a result
			paths->push_back( node->GetPath() );
		}
		else{
			// continue extension
			// all the copies have been used 
			if( newContig->GetOccu() >= newContig->GetCopyNum() )
				continue;
			else{
				newContig->AddOccu();
				string oriString = "BE";
				if( newContigOri == MINUS )
					oriString = "EB";
				char buffer[200];
				sprintf( buffer, "%s\t%s\t%d\t0\n", newContig->GetName().c_str(), oriString.c_str(), newContig->GetLength() );
				string tempPath(buffer);
				PathNode *newNode = new PathNode( newContig, newContigOri, node->GetStep() + 1, node->GetDis() + newContig->GetLength() - Configure::OVERLAP, 
						node->GetPath() + tempPath );

				if( !FindPath( newNode, paths, gap, false ) ){
					newContig->RemoveOccu();
					delete newNode;
					return false;
				}

				newContig->RemoveOccu();
				delete newNode;
			}
		}
	}

	if( paths->empty() && isStart )
	{
		m_disconnectedGap++;
	}

	return true;
}

// select proper path to fill the gap
// find the path which has minimum error from gap estimation
void ContigGraphAnalyzer::SeletePath( Gap *gap, list<string> *paths )
{
	int error = 99999;
	list<string>::iterator bestPath;

	list<string> *lines = new list<string>;
	vector<string> *column = new vector<string>;
	for( list<string>::iterator pathIter = paths->begin(); pathIter != paths->end(); pathIter++ )
	{
		int sum = 0;
		Split( *pathIter, "\n", lines );
		for( list<string>::iterator lineIter = lines->begin(); lineIter != lines->end(); lineIter++ ){
			Split( *lineIter, "\t", column );
			sum += atoi( (*column)[ 2 ].c_str() ) - Configure::OVERLAP;
		}
		if( abs( sum - gap->GetGapSize() ) < error ){
			// find better solution
			bestPath = pathIter;
			error = abs( sum - gap->GetGapSize() );
		}
	}

	// handle the best path, add occurrence for each contig on the path
	Split( *bestPath, "\n", lines );
	vector<string> *names = new vector<string>;
	for( list<string>::iterator lineIter = lines->begin(); lineIter != lines->end(); lineIter++ ){
		Split( *lineIter, "\t", column );
		Split( (*column)[ 0 ], "_", names );
		Contig *c;
		if( m_fileType == VELVET )
			c = (*m_contigs)[ atoi( (*names)[ 1 ].c_str() ) ];
		else
			c = (*m_contigsName)[ (*column)[ 0 ] ];
		if( c != gap->GetStartContig() && c != gap->GetEndContig() )
			c->AddOccu();
	}

	//cout<<*bestPath<<endl;
	if( (*bestPath) == "" )
		gap->SetGap( "Connect" );
	else
		gap->SetGap( *bestPath );

	lines->clear();
	delete lines;
	column->clear();
	delete column;
	names->clear();
	delete names;
}

// read contig fasta file
int ContigGraphAnalyzer::ReadContigFasta( string fileName )
{
	ifstream fastaReader( fileName.c_str() );

	if( fastaReader == NULL ){
		cout<<"\tread contig fasta file error\n";
		return -1;
	}

	string line;
	vector<string> *column = new vector<string>;
	while( getline( fastaReader, line ) ){
		if( line.substr( 0, 1 ) == ">" ){
			// save name
			string name = line.substr( 1, line.length() - 1 );
			Split( line, "_", column );
			int id = atoi( (*column)[ 1 ].c_str() );
			(*m_contigs)[ id ]->SetName( name );
		}
	}

	column->clear();
	delete column;
	fastaReader.close();
	return 1;
}

// analyze subgraph for each gaps
void ContigGraphAnalyzer::AnalyzeGaps( string scafFileName  )
{
	// read scaffold file
	ifstream scafReader( scafFileName.c_str() );

	if( scafReader == NULL )
		return;

	// read all scaffold
	string preline = "", nextline;
	int gapID = 0;
	vector<Gap*> *gaps = new vector<Gap*>();
	list<Gap*> *sortedGaps = new list<Gap*>();

	while( getline( scafReader, nextline ) ){
		if( nextline.substr( 0, 1 ) == ">" ){
			// save last contig
			if( preline != "" )
			{
				vector<string> *temp = new vector<string>;
				Split( preline, "\t", temp );
				Contig *currentContig;
				if( m_fileType == SOAP )
					currentContig = (*m_contigsName)[ (*temp)[ 0 ] ];
				else if( m_fileType == VELVET ){
					vector<string> *tempName = new vector<string>;
					Split( (*temp)[ 0 ], "_", tempName );
					currentContig = (*m_contigs)[ atoi( (*tempName)[ 1 ].c_str() ) ];
					tempName->clear();
					delete tempName;
				}
				else if( m_fileType == ABYSS ){
					currentContig = (*m_contigsName)[ (*temp)[ 0 ] ];
				}
				else{
					cout<<"File Type error!"<<endl;
					return;
				}
				temp->clear();
				delete temp;

				currentContig->AddOccu();
				currentContig->SetInScaffold();
				currentContig->SetCopyNum( 1 );
			}

			// start a new scaffold
			getline( scafReader, preline );
		}
		else{
			vector<string> *temp = new vector<string>;
			Split( preline, "\t", temp );
			Contig *currentContig;
			if( m_fileType == SOAP )
				currentContig = (*m_contigsName)[ (*temp)[ 0 ] ];
			else if( m_fileType == VELVET ){
				vector<string> *tempName = new vector<string>;
				Split( (*temp)[ 0 ], "_", tempName );
				currentContig = (*m_contigs)[ atoi( (*tempName)[ 1 ].c_str() ) ];
				tempName->clear();
				delete tempName;
			}
			else if( m_fileType == ABYSS ){
				currentContig = (*m_contigsName)[ (*temp)[ 0 ] ];
			}
			else{
				cout<<"File Type error!"<<endl;
				return;
			}
			temp->clear();
			delete temp;

			currentContig->AddOccu();
			currentContig->SetInScaffold();
			currentContig->SetCopyNum( 1 );

			// generate a gap
			Gap *newGap = GenerateGap( preline, nextline );
			newGap->SetID( gapID );
			gaps->push_back( newGap );
			sortedGaps->push_back( newGap );
			gapID++;

			preline = nextline;
		}
	}
	scafReader.close();

#ifdef DEBUG
	//cout<<"Total gaps is "<<gapID<<endl;
#endif
	
	// for each gap, get the subgraph 
	for( vector<Gap*>::iterator gapIter = gaps->begin(); gapIter != gaps->end(); gapIter++ )
	{
		GenSubgraph( *gapIter );
	}


	// print connected graph to produce quadratic formulation
	// and solve them using quadratic programming method.
	// Then output all the gaps which cannot be solved.
	//PrintConnectedGaps( gaps );


	////////////Fill the gaps
	///////////Use simple method
	// sort the gaps according to size
	sortedGaps->sort( small_gap );

	// handle each gap in greedy fation accordingly
	ofstream gapWriter( "gapSeq.fasta" );
	for( list<Gap*>::iterator iter = sortedGaps->begin(); iter != sortedGaps->end(); iter++ )
	{
		cout<<(*iter)->GetGapSize()<<endl;
		FillGap( *iter );
		if( (*iter)->GetCalculatedGapSize() != -999 )
		{
			string name = ">" + (*iter)->GetStartString() + "\t" + itos( (*iter)->GetCalculatedGapSize() ) + "\n";
			gapWriter.write( name.c_str(), name.length() );
			gapWriter.write( (*iter)->GetGapSeq().c_str(), (*iter)->GetGapSeq().length() );
		}
	}
	gapWriter.close();

	cout<<"Total gaps: "<<sortedGaps->size()<<endl;
	cout<<"Total filled gaps: "<<this->m_totalFilledGaps<<endl;
	cout<<"# of gaps with multiple paths: "<<this->m_multiPath<<endl;


	vector<Gap*>::iterator iter = gaps->begin();
	while( iter != gaps->end() )
	{
		delete (*iter);
		iter = gaps->erase( iter );
	}
	gaps->clear();
	delete gaps;
	sortedGaps->clear();
	delete sortedGaps;
}

// print connected gaps
void ContigGraphAnalyzer::PrintConnectedGaps( vector<Gap*> *gaps )
{
	int num = 1;
	m_simplePath = "";
	for( vector<Gap*>::iterator gapIter = gaps->begin(); gapIter != gaps->end(); gapIter++ )
	{
		Gap* currentGap = *gapIter;
		int numOfGap = 1;
		int numOfVariables = 0;
		if( !currentGap->IfPrint() )
		{ // print this gap
			set<Gap*> *connectedGaps = new set<Gap*>;
			set<Contig*> *subgraph = new set<Contig*>;
			list<Gap*> *relatedGaps = new list<Gap*>;

			relatedGaps->push_back( currentGap );
			currentGap->SetIfPrint( true );
			connectedGaps->insert( currentGap );

			while( !relatedGaps->empty() )
			{
				Gap *newGap = relatedGaps->front();
				relatedGaps->pop_front();
				numOfVariables += newGap->GetNumOfEdges();
				cout<<newGap->GetID()<<endl;

				subgraph->insert( newGap->GetStartContig() );
				subgraph->insert( newGap->GetEndContig() );

				set<Contig*> *nodes = newGap->GetSubgraph();
				for( set<Contig*>::iterator conIter = nodes->begin(); conIter != nodes->end(); conIter++ )
				{
					subgraph->insert( *conIter );
					list<Gap*>* newGaps = (*conIter)->GetRelatedGaps();
					for( list<Gap*>::iterator gIter = newGaps->begin(); gIter != newGaps->end(); gIter++ )
					{
						if( !(*gIter)->IfPrint() )
						{
							relatedGaps->push_back( *gIter );
							connectedGaps->insert( *gIter );
							(*gIter)->SetIfPrint( true );
							numOfGap++;
						}
					}
				}
			}

			// simplify the graph by combine bubbles and simple paths
			SimplifyGraph( subgraph );

			// record all the simple path
			for( set<Contig*>::iterator nodeIter = subgraph->begin(); nodeIter != subgraph->end(); nodeIter++ )
			{
				if( (*nodeIter)->GetCombinedNodes()->size() > 1 )
				{
					// generate plus path
					m_simplePath += (*nodeIter)->GetName() + "=";
					for( list<pair<string, string> >::iterator nameIter = (*nodeIter)->GetCombinedNodes()->begin(); nameIter != (*nodeIter)->GetCombinedNodes()->end(); nameIter++ )
					{
						if( (*nameIter).first == "-" )
							m_simplePath += (*nameIter).first + (*nameIter).second + "\t";
						else
							m_simplePath += (*nameIter).second + "\t";
					}
					m_simplePath += "\n";

					// generate minus path
					m_simplePath += "-" + (*nodeIter)->GetName() + "=";
					for( list<pair<string, string> >::reverse_iterator nameIter = (*nodeIter)->GetCombinedNodes()->rbegin(); nameIter != (*nodeIter)->GetCombinedNodes()->rend(); nameIter++ )
					{
						if( (*nameIter).first == "+" )
							m_simplePath += "-" + (*nameIter).second + "\t";
						else
							m_simplePath += (*nameIter).second + "\t";
					}
					m_simplePath += "\n";
				}
			}

			int numOfEdge = PrintSubgraph( num, subgraph, "connectedGaps" );

			cout<<"the number of gaps in connectedGaps "<<(num-1)<<"\t"<<numOfGap<<endl;
			cout<<"Number of edges:\t"<<numOfEdge<<endl;
			cout<<"Number of variables:\t"<<numOfVariables<<endl;

			cout<<"generating the quadratic programs..."<<endl;
			if( numOfEdge > 0 ){
				GenQuadProg( connectedGaps, num );
				num++;
			}

			subgraph->clear();
			delete subgraph;
			relatedGaps->clear();
			delete relatedGaps;
			connectedGaps->clear();
			delete connectedGaps;
		}
	}

	OutputSimplePath( "simplePath" );
}

void ContigGraphAnalyzer::OutputSimplePath( string fileName )
{
	ofstream pathWriter( fileName.c_str() );

	if( pathWriter == NULL ){
		cout<<"error writing "<<fileName<<endl;
		return;
	}

	pathWriter.write( m_simplePath.c_str(), m_simplePath.length() );
}

// simplify the connected graph
void ContigGraphAnalyzer::SimplifyGraph( set<Contig*> *nodes )
{
	bool handleBubble;
	bool handleSimplePath;
	do{
		handleBubble = false;
		handleSimplePath = false;
		
		// check if there is a bubble
		bool hasBubble;
		do{
			hasBubble = CheckBubble( nodes );
			if( hasBubble ){
				handleBubble = true;
			}
		}while( hasBubble );

		// check if there is a simple path
		bool hasSimplePath;
		do{
			hasSimplePath = CheckSimplePath( nodes );
			if( hasSimplePath )
				handleSimplePath = true;

		}while( hasSimplePath );
	}while( handleBubble || handleSimplePath );
}

// check if the graph contain simple path
bool ContigGraphAnalyzer::CheckSimplePath( set<Contig*> *nodes )
{
	for( set<Contig*>::iterator nodeIter = nodes->begin(); nodeIter != nodes->end(); nodeIter++ )
	{
		Contig *c = *nodeIter;
		// do not combine scaffold node
		if( c->IfInScaffold() )
			continue;
		
		if( IfBelongToSimplePath( c,  nodes, LEFT ) ){
			printf( "find a simple path: %s left\n", c->GetName().c_str() );
			return true;
		}

		if( IfBelongToSimplePath( c,  nodes, RIGHT ) ){
			printf( "find a simple path: %s right\n", c->GetName().c_str() );
			return true;
		}
	}
	return false;
}

// check if current contig c belong to a simple path
bool ContigGraphAnalyzer::IfBelongToSimplePath( Contig *c, set<Contig*> *nodes, int ori)
{
	// check if belong to simple path
	Contig *otherContig;
	int numOfEdgeC = 0;
	int numOfEdgeOther = 0;
	int extentOriOther;			// the extention orientation of other contig on the path
	int keyOfOtherInC;

	// check if c has only one edge
	map<int, pair<Contig*, int> > *edges = c->GetEdge( ori );
	for( map<int, pair<Contig*, int> >::iterator edgeIter = edges->begin(); edgeIter != edges->end(); edgeIter++ )
	{
		if( nodes->find( (*edgeIter).second.first ) != nodes->end() ){
			numOfEdgeC++;
			if( numOfEdgeC > 1 )
				break;
			
			otherContig = (*edgeIter).second.first;
			keyOfOtherInC = (*edgeIter).first;

			if( (*edgeIter).first > 0 ){
				if( ori == LEFT )
					extentOriOther = RIGHT;
				else
					extentOriOther = LEFT;
			}
			else{
				if( ori == LEFT )
					extentOriOther = LEFT;
				else
					extentOriOther = RIGHT;
			}

		}
	}
	if( numOfEdgeC == 1 )
	{ // find a contig
		// do not consider scaffold node
		if( otherContig->IfInScaffold() )
			return false;

		// check if other contig also only has one connection
		edges = otherContig->GetEdge( extentOriOther );
		int keyOfCInOther;
		pair<Contig*, int> valueOfCInOther;
		int num = 0;
		for( map<int, pair<Contig*, int> >::iterator edgeIter = edges->begin(); edgeIter != edges->end(); edgeIter++ )
		{
			if( nodes->find( (*edgeIter).second.first ) != nodes->end() ){
				num++;
				keyOfCInOther = (*edgeIter).first;
				valueOfCInOther = (*edgeIter).second;
			}
		}

		if( num == 1 )
		{ // it is a simple path, combine
			// update the combined nodes into c
			if( ori == RIGHT )
			{ // extend to right
				if( extentOriOther == LEFT )
				{ // + +
					c->AddCombinedNode( otherContig, ori, PLUS );
				}
				else
				{ // + -
					c->AddCombinedNode( otherContig, ori, MINUS );
				}
			}
			else
			{ // extend to left
				if( extentOriOther == LEFT )
				{ // - +
					c->AddCombinedNode( otherContig, ori, MINUS );
				}
				else
				{ // - -
					c->AddCombinedNode( otherContig, ori, PLUS );
				}
			}

			// remove other node
			printf( "remove %s\n", otherContig->GetName().c_str() );
			nodes->erase( otherContig );
			otherContig->RemoveFromAllGaps();

			// update the length of new contig 
			c->SetLength( c->GetLength() + otherContig->GetLength() - Configure::OVERLAP );

			// add the edges of other node to contig c
			int oriOfOtherToMove;		// the direction of extention of other contig to move the edges
			
			if( extentOriOther == LEFT ){
				printf( "extentOriOther is LEFT\n" );
				oriOfOtherToMove = RIGHT;
			}
			else{
				printf( "extentOriOther is RIGHT\n" );
				oriOfOtherToMove = LEFT;
			}

			edges = otherContig->GetEdge( oriOfOtherToMove );
			for( map<int, pair<Contig*, int> >::iterator edgeIter = edges->begin(); edgeIter != edges->end(); edgeIter++ )
			{
				printf( "moving edge of %s\n", (*edgeIter).second.first->GetName().c_str() );
				if( ori == oriOfOtherToMove ){
					printf( "same ori\n" );
					printf( "%d\n", (*edgeIter).first );
				}
				if( ori == oriOfOtherToMove )
					(*c->GetEdge( ori ))[ (*edgeIter).first ] = (*edgeIter).second;
				else
					(*c->GetEdge( ori ))[ -(*edgeIter).first ] = (*edgeIter).second;

				int newOri;
				if( oriOfOtherToMove == RIGHT ){
					if( (*edgeIter).first > 0 )
						newOri = LEFT;
					else
						newOri = RIGHT;
				}
				else{
					if( (*edgeIter).first > 0 )
						newOri = RIGHT;
					else
						newOri = LEFT;
				}

				if( newOri == LEFT )
					printf( "newOri: LEFT\n" );
				else
					printf( "newOri: Right\n" );

				if( (*edgeIter).first * keyOfCInOther > 0 )
					printf( "same key sign\n" );
				
				printf( "%d\n", keyOfCInOther );
				if( ori == oriOfOtherToMove ){ 
					if( (*edgeIter).first * keyOfCInOther > 0 )
						(*(*edgeIter).second.first->GetEdge( newOri ))[ keyOfCInOther ] = valueOfCInOther;
					else
						(*(*edgeIter).second.first->GetEdge( newOri ))[ -keyOfCInOther ] = valueOfCInOther;
				}
				else
					if( (*edgeIter).first * keyOfCInOther > 0 )
						(*(*edgeIter).second.first->GetEdge( newOri ))[ -keyOfCInOther ] = valueOfCInOther;
					else
						(*(*edgeIter).second.first->GetEdge( newOri ))[ keyOfCInOther ] = valueOfCInOther;
			}
			return true;
		}
	}

	return false;
}

// check if the graph contain bubble
// find the first bubble, combine it and return;
bool ContigGraphAnalyzer::CheckBubble( set<Contig*> *nodes )
{
	for( set<Contig*>::iterator nodeIter = nodes->begin(); nodeIter != nodes->end(); nodeIter++ )
	{
		Contig *currentContig = *nodeIter;
		if( currentContig->IfInScaffold() ) 
			continue;

		if( IfBelongToBubble( currentContig,  nodes ) ){
			//printf( "find a bubble: %s\n", currentContig->GetName().c_str() );
			return true;
		}
	}
	return false;
}

// check if current contig c belong to a bubble
bool ContigGraphAnalyzer::IfBelongToBubble( Contig *c, set<Contig*> *nodes )
{
	int leftSize = 0;
	int rightSize = 0;
	int leftCKey = 0;
	int rightCKey = 0;
	int leftOtherKey = 0;
	int rightOtherKey = 0;
	int numOfLeftEdgeC = 0;
	int numOfRightEdgeC = 0;
	int cClusterSize = 0;		// the cluster size of c;
	int rightOri;
	int leftOri;
	Contig *leftContig, *rightContig;

	int startKeyOfC, endKeyOfC;		// the key of contig c in start and end contig

	// check if c has only one out edge and only one in edge
	map<int, pair<Contig*, int> > *edges = c->GetEdge( LEFT );
	for( map<int, pair<Contig*, int> >::iterator edgeIter = edges->begin(); edgeIter != edges->end(); edgeIter++ )
	{
		if( nodes->find( (*edgeIter).second.first ) != nodes->end() ){
			numOfLeftEdgeC++;
			if( numOfLeftEdgeC > 1 )
				return false;
			cClusterSize += (*edgeIter).second.second;
			leftSize = (*edgeIter).second.second;
			leftCKey = (*edgeIter).first;
			leftContig = (*edgeIter).second.first;
			if( (*edgeIter).first > 0 )
				leftOri = RIGHT;
			else
				leftOri = LEFT;
		}
	}
	if( numOfLeftEdgeC == 0 )
		return false;

	edges = c->GetEdge( RIGHT );
	for( map<int, pair<Contig*, int> >::iterator edgeIter = edges->begin(); edgeIter != edges->end(); edgeIter++ )
	{
		if( nodes->find( (*edgeIter).second.first ) != nodes->end() ){
			numOfRightEdgeC++;
			if( numOfRightEdgeC > 1 )
				return false;
			cClusterSize += (*edgeIter).second.second;
			rightSize = (*edgeIter).second.second;
			rightCKey = (*edgeIter).first;
			rightContig = (*edgeIter).second.first;
			if( edges->begin()->first > 0 )
				rightOri = LEFT;
			else
				rightOri = RIGHT;
		}
	}
	if( numOfRightEdgeC == 0 )
		return false;

 	// check if the start and end contig contain another contig c' like c
	map<int, pair<Contig*, int> > *startNeighbor = leftContig->GetEdge( leftOri );
	map<int, pair<Contig*, int> > *endNeighbor = rightContig->GetEdge( rightOri );

	// get key of c
	for( map<int, pair<Contig*, int> >::iterator neighborIter = startNeighbor->begin(); neighborIter != startNeighbor->end(); neighborIter++ )
	{
		if( (*neighborIter).second.first == c ){
			startKeyOfC = (*neighborIter).first;
			break;
		}
	}
	for( map<int, pair<Contig*, int> >::iterator neighborIter = endNeighbor->begin(); neighborIter != endNeighbor->end(); neighborIter++ )
	{
		if( (*neighborIter).second.first == c ){
			endKeyOfC = (*neighborIter).first;
			break;
		}
	}

	for( map<int, pair<Contig*, int> >::iterator neighborIter = startNeighbor->begin(); neighborIter != startNeighbor->end(); neighborIter++ )
	{
		int nodeID = (*neighborIter).first;
		int leftOtherKey = (*neighborIter).first;
		int startKeyInOther, endKeyInOther, startOriInOther, endOriInOther;
		if( leftOri == rightOri )
			nodeID = -nodeID;

		if( endNeighbor->find( nodeID ) != endNeighbor->end() )
		{ // find a node
			// if c' has similar length as c, combine them
			Contig *otherNode = (*neighborIter).second.first;

			// do not shrink node in scaffold
			if( otherNode->IfInScaffold() )
				continue;
			if( otherNode == c )
				continue;
			if( nodes->find( otherNode ) == nodes->end() )
				continue;
			if( abs( otherNode->GetLength() - c->GetLength() ) / c->GetLength() >= 0.1 )
				continue;

			// check if c has only one out edge and only one in edge
			int numOfLeftEdgeOther = 0;
			int numOfRightEdgeOther = 0;
			int otherClusterSize = 0;

			edges = otherNode->GetEdge( LEFT );
			for( map<int, pair<Contig*, int> >::iterator edgeIter = edges->begin(); edgeIter != edges->end(); edgeIter++ )
			{
				if( nodes->find( (*edgeIter).second.first ) != nodes->end() ){
					if( (*edgeIter).second.first == leftContig ){
						startKeyInOther = (*edgeIter).first;
						startOriInOther = LEFT;
					}
					else if( (*edgeIter).second.first == rightContig ){
						endKeyInOther = (*edgeIter).first;
						endOriInOther = LEFT;
					}
					numOfLeftEdgeOther++;
					if( numOfLeftEdgeOther > 1 )
						return false;
					otherClusterSize += (*edgeIter).second.second;
					leftSize += (*edgeIter).second.second;
				}
			}
			if( numOfLeftEdgeOther == 0 )
				return false;

			edges = otherNode->GetEdge( RIGHT );
			for( map<int, pair<Contig*, int> >::iterator edgeIter = edges->begin(); edgeIter != edges->end(); edgeIter++ )
			{
				if( nodes->find( (*edgeIter).second.first ) != nodes->end() ){
					if( (*edgeIter).second.first == leftContig ){
						startKeyInOther = (*edgeIter).first;
						startOriInOther = RIGHT;
					}
					else if( (*edgeIter).second.first == rightContig ){
						endKeyInOther = (*edgeIter).first;
						endOriInOther = RIGHT;
					}
					numOfRightEdgeOther++;
					if( numOfRightEdgeOther > 1 )
						return false;
					otherClusterSize += (*edgeIter).second.second;
					rightSize += (*edgeIter).second.second;
				}
			}
			if( numOfRightEdgeOther == 0 )
				return false;

			int rightOtherKey = endNeighbor->find( nodeID )->first;
			//printf( "other node is %s\n", otherNode->GetName().c_str() );
			//printf( "leftSize is %d\n", leftSize );
			//printf( "other size is %d, while c size is %d\n", otherClusterSize, cClusterSize );

			if( otherClusterSize > cClusterSize )
			{ // discard node c
				printf( "discard c\n" );
				nodes->erase( c );	// remove from graph

				// update the copy number
				otherNode->SetCopyNum( otherNode->GetCopyNum() + c->GetCopyNum() );

				// remove from gap
				c->RemoveFromAllGaps();
				// update the size
				(*otherNode->GetEdge( startOriInOther ))[ startKeyInOther ] = pair<Contig*, int>(leftContig, leftSize);
				(*otherNode->GetEdge( endOriInOther ))[ endKeyInOther ] = pair<Contig*, int>(rightContig, rightSize);
				(*leftContig->GetEdge( leftOri ))[ leftOtherKey ] = pair<Contig*, int>(otherNode, leftSize);
				(*rightContig->GetEdge( rightOri ))[ rightOtherKey ] = pair<Contig*, int>(otherNode, rightSize);
				return true;
			}
			else
			{ // discard other node
				printf( "discard other node\n" );
				int num = nodes->erase( otherNode );
				printf( "erased %d nodes\n", num );

				// update the copy number
				c->SetCopyNum( c->GetCopyNum() + otherNode->GetCopyNum() );

				otherNode->RemoveFromAllGaps();
				(*c->GetEdge( LEFT ))[ leftCKey ] = pair<Contig*, int>(leftContig, leftSize);
				(*c->GetEdge( RIGHT ))[ rightCKey ] = pair<Contig*, int>(rightContig, rightSize);
				(*leftContig->GetEdge( leftOri ))[ startKeyOfC ] = pair<Contig*, int>(c, leftSize);
				(*rightContig->GetEdge( rightOri ))[ endKeyOfC ] = pair<Contig*, int>(c, rightSize);
				
				return true;
			}
		}
	}

	return false;
}

// generate the subgraph for each gap
void ContigGraphAnalyzer::GenSubgraph( Gap *gap )
{
	cout<<"Gap size is "<<gap->GetGapSize()<<endl;

	// get all connected nodes within maximum distance from the start contig
	set<ConnectedNode*> *startSet = new set<ConnectedNode*>;
	GetPossibleNodes( gap, START, startSet );
	// get all connected nodes within maximum distance from the end contig
	set<ConnectedNode*> *endSet = new set<ConnectedNode*>;
	GetPossibleNodes( gap, END, endSet );

	/*cout<<"current gap is:"<<endl;
	cout<<gap->GetStartString()<<endl<<gap->GetEndString()<<endl;
	cout<<"There are "<<startSet->size()<<" contigs in start set\n";
	for( set<ConnectedNode*>::iterator setIter = startSet->begin(); setIter != startSet->end(); setIter++ )
	{
		cout<<(*setIter)->GetContig()->GetName()<<"\t"<<(*setIter)->GetContig()->GetCov()<<"\tori\t"<<(*setIter)->GetOri()<<"\tdis:"<<(*setIter)->GetMinDis()<<endl;
	}
	cout<<"There are "<<endSet->size()<<" contigs in end set\n";
	for( set<ConnectedNode*>::iterator setIter = endSet->begin(); setIter != endSet->end(); setIter++ )
	{
		cout<<(*setIter)->GetContig()->GetName()<<"\t"<<(*setIter)->GetContig()->GetCov()<<"\tori\t"<<(*setIter)->GetOri()<<"\tdis:"<<(*setIter)->GetMinDis()<<endl;
	}*/

	set<Contig*> *subgraph = new set<Contig*>;		// all the possible nodes in this gap

	// find the possible nodes, L1 + L2 < maximum distance
	map<pair<Contig*, int>, ConnectedNode*> *endMap = new map<pair<Contig*, int>, ConnectedNode*>;
	for( set<ConnectedNode*>::iterator setIter = endSet->begin(); setIter != endSet->end(); setIter++ )
	{
		(*endMap)[ pair<Contig*, int>( (*setIter)->GetContig(), (*setIter)->GetOri() ) ] = *setIter;
	}

	for( set<ConnectedNode*>::iterator setIter = startSet->begin(); setIter != startSet->end(); setIter++ )
	{
		map<pair<Contig*, int>, ConnectedNode*>::iterator mapIter = endMap->find( pair<Contig*, int>( (*setIter)->GetContig(), (*setIter)->GetOri() ) );
		if( mapIter != endMap->end() )
		{ // this node is in the intersection set
			if( (*mapIter).second->GetMinDis() + (*setIter)->GetMinDis() - (*setIter)->GetContig()->GetLength() < gap->GetMaxGapSize() ){
				// sum of two nearest distance is smaller than the upperbound, add to the subgraph
				subgraph->insert( (*setIter)->GetContig() );
			}
		}
	}

	// set the connection 
	cout<<"number of contigs in subgraph: "<<subgraph->size()<<endl;
	for( set<Contig*>::iterator setIter = subgraph->begin(); setIter != subgraph->end(); setIter++ )
	{
		(*setIter)->AddGap( gap );
	}
	gap->SetSubgraph( subgraph );

	// print subgraph
	subgraph->insert( gap->GetStartContig() );
	subgraph->insert( gap->GetEndContig() );
	int numOfEdge = this->PrintSubgraph( gap->GetID(), subgraph, "gap" );
	gap->SetNumOfEdges( numOfEdge );

	// delete all variables
	endMap->clear();
	delete endMap;
	subgraph->clear();
	delete subgraph;
	for( set<ConnectedNode*>::iterator setIter = endSet->begin(); setIter != endSet->end(); setIter++ )
	{
		delete *setIter;
	}
	endSet->clear();
	delete endSet;

	for( set<ConnectedNode*>::iterator setIter = startSet->begin(); setIter != startSet->end(); setIter++ )
	{
		delete *setIter;
	}
	startSet->clear();
	delete startSet;
}


// search the possible nodes in iterative way
// use dijstra algorithm
void ContigGraphAnalyzer::GetPossibleNodes( Gap *gap, int startContigType, set<ConnectedNode*> *nodeSet )
{
	Contig *startContig;
	int startOri;

	if( startContigType == START ){
		startContig = gap->GetStartContig();
		startOri = gap->GetStartOri();
	}
	else{
		startContig = gap->GetEndContig();
		startOri = gap->GetEndOri();
	}

	int maxDis = gap->GetMaxGapSize();

	// Simple algorithm to calculate the smallest distance of all node from start contig
	// S = empty
	// Q = start contig
	// while ( Q is not empty )
	// {
	//		u = smallest in Q;
	//		S = S+{u}
	//		for each vertex of v, which is the neighbour of u and has not been visited
	//			Relax v; Q = Q + {v}
	// }
	vector<ConnectedNode*> *nodeHeap = new vector<ConnectedNode*>;
	ConnectedNode *newNode = new ConnectedNode( startContig, startOri, 0 );
	nodeHeap->push_back( newNode );
	// create minimum heap
	make_heap( nodeHeap->begin(), nodeHeap->end(), smallerDistance() );

	set<pair<Contig*, int> > *visitedNodes = new set<pair<Contig*, int> >;
	bool isFirst = true;		// label the first contig
	while( !nodeHeap->empty() )
	{ 
		// get smallest node
		ConnectedNode *smallestNode = nodeHeap->front();
		pop_heap( nodeHeap->begin(), nodeHeap->end(), smallerDistance() );
		nodeHeap->pop_back();

		Contig *currentContig = smallestNode->GetContig();
		// do not add first contig into the set
		if( isFirst ){
			isFirst = false;
		}
		else{
			nodeSet->insert( smallestNode );
		}

		//visitedNodes->insert( pair<Contig*, int> (currentContig, smallestNode->GetOri() ) );

		// add the neighbours
		map<int, pair<Contig*, int> > *edges;
		int edgeOri;
		if( startContigType == START ){
			// extend to right direction
			if( smallestNode->GetOri() == PLUS )
				edgeOri = RIGHT;
			else
				edgeOri = LEFT;
		}
		else{
			// extend to left direction
			if( smallestNode->GetOri() == PLUS )
				edgeOri = LEFT;
			else
				edgeOri = RIGHT;
		}
		edges = currentContig->GetEdge( edgeOri );

		for( map<int, pair<Contig*, int> >::iterator edgeIter = edges->begin(); edgeIter != edges->end(); edgeIter++ )
		{
			Contig *newContig = (*edgeIter).second.first;
			int newOri = (*edgeIter).first;
			
			int overlap = Configure::OVERLAP;
			if( m_fileType == ABYSS )
				overlap = (*edgeIter).second.second;				// for abyss, this position saves the overlap length

			// do not consider contigs with low coverage
			if( newContig->GetCov() < m_covThreshold )
				continue;

			// if current contig has minus orientation, the following orientation of new contig has to be reversed
			if( smallestNode->GetOri() == PLUS ){
				if( newOri > 0 )
					newOri = PLUS;
				else
					newOri = MINUS;
			}
			else{
				if( newOri > 0 )
					newOri = MINUS;
				else
					newOri = PLUS;
			}

			int maxFutureOverlap = Configure::OVERLAP;
			if( m_fileType == ABYSS ){
				// get the biggest overlap of future extension
				map<int, pair<Contig*, int> > *futureEdges;
				int edgeOri;
				if( startContigType == START ){
					// extend to right direction
					if( newOri == PLUS )
						edgeOri = RIGHT;
					else
						edgeOri = LEFT;
				}	
				else{
					// extend to left direction
					if( newOri == PLUS )
						edgeOri = LEFT;
					else
						edgeOri = RIGHT;
				}
				futureEdges = newContig->GetEdge( edgeOri );

				maxFutureOverlap = 0;
				for( map<int, pair<Contig*, int> >::iterator futureEdgeIter = futureEdges->begin(); futureEdgeIter != futureEdges->end(); futureEdgeIter++ )
				{
					if( (*futureEdgeIter).second.second > maxFutureOverlap )
						maxFutureOverlap = (*futureEdgeIter).second.second;
				}
			}

			// do not consider the contigs which occurs more than copy number
			if( newContig->GetOccu() >= newContig->GetCopyNum() )
				continue;

			// do not consider too long path
			if( smallestNode->GetMinDis() + newContig->GetLength() - overlap - maxFutureOverlap > maxDis )
				continue;

			if( visitedNodes->find( pair<Contig*, int> (newContig, newOri) ) == visitedNodes->end() )
			{ // if this node has not been visited, add to set
				ConnectedNode *newConnectedNode = new ConnectedNode( newContig, newOri, smallestNode->GetMinDis() + newContig->GetLength() - overlap );
				nodeHeap->push_back( newConnectedNode );
				push_heap( nodeHeap->begin(), nodeHeap->end(), smallerDistance() );

				visitedNodes->insert( pair<Contig*, int> (newContig, newOri ) );
			}
			else
			{ // relax

			}
		}
	}

	nodeHeap->clear();
	delete nodeHeap;
	visitedNodes->clear();
	delete visitedNodes;
}

// read SOAPdenovo file to get kmer size
int ContigGraphAnalyzer::GetKmerOfSOAP( string prefix )
{
	string kmerFile = prefix + ".preGraphBasic";

	// read scaffold file
	ifstream scafReader( kmerFile.c_str() );

	if( scafReader == NULL ){
		cout<<"Error in reading \""<<kmerFile<<"\""<<endl;
		return -1;
	}

	// read all scaffold
	string line, nextline;
	getline( scafReader, line );
	vector<string> *temp = new vector<string>;
	Split( line, " ", temp );
	Configure::OVERLAP = atoi( (*temp)[ 3 ].c_str() );
	//cout<<"kmer: "<<m_kmer<<endl;
	temp->clear();
	delete temp;

	return 0;
}

// print graph around vector
int ContigGraphAnalyzer::PrintGraphAroundVector()
{
	set<Contig*> *vectors = new set<Contig*>;

	// read mapping file
	cout<<"reading mapping file"<<endl;
	string mappingFileName = "vector.coords";
	vector<string> *line = new vector<string>;

	ifstream mapReader( mappingFileName.c_str() );
	if( mapReader == NULL )
	{
		cout<<"Error reading mapping file: "<<mappingFileName<<endl;
		delete line;
		delete vectors;
		return -1;
	}

	// read mapping file to recognize the vector contigs
	string lineStr;
	for( int i = 0; i < 4; i++ )
		getline( mapReader, lineStr );

	while( getline( mapReader, lineStr ) )
	{
		Split( lineStr, "\t", line );
		Contig *currentContig = (*m_contigsName)[ line->at( 12 ) ];
		vectors->insert( currentContig );

		int startPos = atoi( line->at( 0 ).c_str() );
		int endPos = atoi( line->at( 1 ).c_str() );
		int ori = PLUS;
		if( atoi( line->at( 2 ).c_str() ) > atoi( line->at( 3 ).c_str() ) )
			ori = MINUS;
		currentContig->SetVector( startPos, endPos, ori );
		currentContig->SetIfPrintInVec( false );
		currentContig->SetIfVec( true );
	}
	delete line;

	// print graph around vector contigs
	/*cout<<"analyzing graph file"<<endl;
	int step = 2;		// the threshold of step outside the vector contigs
	int vecID = 1;
	for( set<Contig*>::iterator contigIter = vectors->begin(); contigIter != vectors->end(); contigIter++ )
	{
		Contig *currentContig = *contigIter;
		set<Contig*> *subgraph = new set<Contig*>;
		
		if( !currentContig->IfPrintInVec() ){
			// establish the graph
			currentContig->SetStepInVecGraph( 0 );
			currentContig->SetIfPrintInVec( true );
			subgraph->insert( currentContig );

			list<Contig*> *possibleContigs = new list<Contig*>;
			possibleContigs->push_back( currentContig );

			while( !possibleContigs->empty() )
			{
				Contig *firstContig = possibleContigs->front();
				possibleContigs->pop_front();

				if( firstContig->GetStepInVecGraph() >= step )
					continue;

				map<int, pair<Contig*, int> > *edges = new map<int, pair<Contig*, int> >;
				edges->insert( firstContig->GetEdge( LEFT )->begin(), firstContig->GetEdge( LEFT )->end() );
				edges->insert( firstContig->GetEdge( RIGHT )->begin(), firstContig->GetEdge( RIGHT )->end() );

				for( map<int, pair<Contig*, int> >::iterator edgeIter = edges->begin(); edgeIter != edges->end(); edgeIter++ )
				{
					Contig *nextContig = (*edgeIter).second.first;
					if( !nextContig->IfPrintInVec() )
					{
						nextContig->SetIfPrintInVec( true );
						subgraph->insert( nextContig );
						possibleContigs->push_back( nextContig );

						if( nextContig->IfVec() )
							nextContig->SetStepInVecGraph( 0 );
						else
							nextContig->SetStepInVecGraph( firstContig->GetStepInVecGraph() + 1 );
					}
				}
				edges->clear();
				delete edges;
			}
			
			possibleContigs->clear();
			delete possibleContigs;
		}

		if( subgraph->size() > 1 )
			PrintSubgraph( vecID++, subgraph, "vec" );
		
		subgraph->clear();
		delete subgraph;
	}*/

	// print complete vector graph with 1 step extension
	cout<<"analyzing graph file"<<endl;
	set<Contig*> *subgraph = new set<Contig*>;
	for( set<Contig*>::iterator contigIter = vectors->begin(); contigIter != vectors->end(); contigIter++ )
	{
		Contig *currentContig = *contigIter;

		// add this contig into subgraph
		subgraph->insert( currentContig );
		
		// add all neighbours
		map<int, pair<Contig*, int> > *edges = new map<int, pair<Contig*, int> >;
		edges->insert( currentContig->GetEdge( LEFT )->begin(), currentContig->GetEdge( LEFT )->end() );
		edges->insert( currentContig->GetEdge( RIGHT )->begin(), currentContig->GetEdge( RIGHT )->end() );
		for( map<int, pair<Contig*, int> >::iterator edgeIter = edges->begin(); edgeIter != edges->end(); edgeIter++ )
		{
			Contig *nextContig = (*edgeIter).second.first;
			subgraph->insert( nextContig );
		}
		edges->clear();
		delete edges;
	}

	PrintSubgraph( 1, subgraph, "allVector" );
		
	subgraph->clear();
	delete subgraph;

	vectors->clear();
	delete vectors;

	return 0;
}



void ContigGraphAnalyzer::GenQuadProg( set<Gap*> *gaps, int id ){
	string fileName = m_outputPath + "quadProg" + itos( id ) + ".lp";
	FILE *progFile = fopen( fileName.c_str(), "w" );
	fprintf( progFile, "\\quadProg%d\n", id );

	set<Contig*> *wholegraph = new set<Contig*>;
	vector<int> *variables = new vector<int>;
	variables->push_back( 0 );
	int numOfVar = 0;
	string problemStr = "";
	m_aimStringQuad = "";
	m_aimStringLinear = "";
	m_constraints = "";
	m_variableString = "  ";
	m_aimLinearValue->clear();

	// save the contigs in scaffold
	vector<int> *scafVar = new vector<int>;		// the variables related with contigs in scaffold (start and end)
	for( vector<vector<int>*>::iterator quadIter = m_aimQuadValue->begin(); quadIter != m_aimQuadValue->end(); quadIter++ )
	{
		(*quadIter)->clear();
		delete *quadIter;
	}
	m_aimQuadValue->clear();

	for( set<Gap*>::iterator gapIter = gaps->begin(); gapIter != gaps->end(); gapIter++ )
	{
		scafVar->clear();
		problemStr += "\\One gap graph:\n";
		Gap *currentGap = *gapIter;
		set<Contig*> *graph = new set<Contig*>;
		graph->insert( currentGap->GetSubgraph()->begin(), currentGap->GetSubgraph()->end() );
		graph->insert( currentGap->GetStartContig() );
		graph->insert( currentGap->GetEndContig() );
		
		// handle each edge and create variables
		for( set<Contig*>::iterator contigIter = graph->begin(); contigIter != graph->end(); contigIter++ ){
			Contig *currentContig = (*contigIter);
			wholegraph->insert( currentContig );
			CreateVariable( currentContig, currentContig->GetEdge( LEFT ), LEFT, graph, variables, problemStr, scafVar );
			CreateVariable( currentContig, currentContig->GetEdge( RIGHT ), RIGHT, graph, variables, problemStr, scafVar );
		}

		// add something to aim function
		GenAimFunction( currentGap, problemStr );

		// add punishment to aim function
		//AddPunishFunction( scafVar, variables );

		// construct the constraints about in==out
		// and combine new variables to all variables 
		problemStr += "\\Constraints:\n";
		for( set<Contig*>::iterator contigIter = graph->begin(); contigIter != graph->end(); contigIter++ ){
			Contig *currentContig = (*contigIter);
			if( currentContig != currentGap->GetStartContig() && currentContig != currentGap->GetEndContig() )
				GenInOutEquation( currentContig, problemStr );
			currentContig->CombineVariables();
		}
	}

	// calculate quadratic coefficients
	for( int i = 0; i < m_aimQuadValue->size(); i++ )
	{
		vector<int> *iVector = m_aimQuadValue->at( i );
		if( iVector->size() > i )
			(*iVector)[ i ] = (*iVector)[i] * (*iVector)[i] * 2;
	}

	// construct the constraints about copy number and clear the variables
	problemStr += "\\Copy number constraints:\n";
	for( set<Contig*>::iterator contigIter = wholegraph->begin(); contigIter != wholegraph->end(); contigIter++ )
	{
		if( (*contigIter)->IfInScaffold() )
		{ // it is in scaffold, have to check two directions
			if( !(*contigIter)->GetInEdgeVar()->empty() ){
				if( (*contigIter)->CheckIfExtendLeft(gaps) )
					GenEqualOneConstraints( (*contigIter)->GetInEdgeVar(), *contigIter, problemStr );
				else
					GenCopyNumberConstraints( (*contigIter)->GetInEdgeVar(), *contigIter, problemStr );
				
			}

			if( !(*contigIter)->GetOutEdgeVar()->empty() ){
				if( (*contigIter)->CheckIfExtendRight(gaps) )
					GenEqualOneConstraints( (*contigIter)->GetOutEdgeVar(), *contigIter, problemStr );
				else
					GenCopyNumberConstraints( (*contigIter)->GetOutEdgeVar(), *contigIter, problemStr );
				
			}
		}
		else if( !(*contigIter)->GetInEdgeVar()->empty() )
			GenCopyNumberConstraints( (*contigIter)->GetInEdgeVar(), *contigIter, problemStr );

		(*contigIter)->ClearVar();
	}

	fprintf( progFile, "%s", problemStr.c_str() );
	
	// print aim function
	fprintf( progFile, "min\n" );
	//fprintf( progFile, "myobj: %s + [ %s ] / 2\n", m_aimStringLinear.c_str(), m_aimStringQuad.c_str() );
	fprintf( progFile, "myobj: " );
	if( m_aimLinearValue->size() > 2 ) 
	{
		fprintf( progFile, "%d X1 ", m_aimLinearValue->at( 1 ) );
	}

	for( int i = 2; i < m_aimLinearValue->size(); i++ )
	{
		if( m_aimLinearValue->at( i ) < 0 )
			fprintf( progFile, " - %d X%d", abs(m_aimLinearValue->at( i )), i );
		else
			fprintf( progFile, " + %d X%d", abs(m_aimLinearValue->at( i )), i );
	}

	// print quadratic format
	fprintf( progFile, " + [ " );
	bool firstQuad = true;
	for( int i = 1; i < m_aimQuadValue->size(); i++ )
	{
		vector<int> *row = (*m_aimQuadValue)[ i ];
		for( int j = 1; j < row->size(); j++ )
		{
			if( row->at( j ) != 0  )
			{
				if( firstQuad )
				{
					fprintf( progFile, "%d ", row->at( j ) );
					firstQuad = false;
				}
				else
				{
					if( row->at( j ) > 0 )
						fprintf( progFile, "+ %d ", row->at( j ) );
					else
						fprintf( progFile, "- %d ", abs(row->at( j )) );
				}

				if( i == j )
					fprintf( progFile, "X%d ^ 2 ", i );
				else
					fprintf( progFile, "X%d * X%d ", i, j );
			}
		}
	}
	fprintf( progFile, " ] / 2\n" );

	//fprintf( progFile, " + [ %s ] / 2\n", m_aimStringQuad.c_str() );

	// print constraints
	fprintf( progFile, "st\n%s", m_constraints.c_str() );
	fprintf( progFile, "general\n%s\n", m_variableString.c_str() );
	fprintf( progFile, "end\n" );

	wholegraph->clear();
	delete wholegraph;
	variables->clear();
	delete variables;
	fclose( progFile );
}

void ContigGraphAnalyzer::GenCopyNumberConstraints( list<int> *var, Contig *c, string &problemStr )
{
	list<int>::iterator varIter = var->begin();
	problemStr += "\\X" + itos( *varIter );
	m_constraints += "X" + itos( *varIter );
	varIter++;
	for( ; varIter != var->end(); varIter++ )
	{
		problemStr += "\t+\tX" + itos( *varIter );
		m_constraints += "\t+\tX" + itos( *varIter );
	}

	problemStr += "\t<=\t" + itos( c->GetCopyNum() ) + "\n";
	m_constraints += "\t<=\t" + itos( c->GetCopyNum() ) + "\n";
}

void ContigGraphAnalyzer::GenEqualOneConstraints( list<int> *var, Contig *c, string &problemStr )
{
	list<int>::iterator varIter = var->begin();
	problemStr += "\\X" + itos( *varIter );
	m_constraints += "X" + itos( *varIter );
	varIter++;
	for( ; varIter != var->end(); varIter++ )
	{
		problemStr += "\t+\tX" + itos( *varIter );
		m_constraints += "\t+\tX" + itos( *varIter );
	}

	problemStr += "\t=\t1\n";
	m_constraints += "\t=\t1\n";
}

void ContigGraphAnalyzer::GenInOutEquation( Contig *c, string &problemStr )
{
	list<int> *edges = c->GetNewInEdgeVar();
	list<int>::iterator varIter = edges->begin();
	problemStr += "\\X" + itos( *varIter );
	m_constraints += "X" + itos( *varIter );
	varIter++;
	while( varIter != edges->end() )
	{
		problemStr += "\t+\tX" + itos( *varIter );
		m_constraints += "\t+\tX" + itos( *varIter );
		varIter++;
	}

	problemStr += "\t=\t";

	edges = c->GetNewOutEdgeVar();
	varIter = edges->begin();
	problemStr += "X" + itos( *varIter );
	m_constraints += "\t-\tX" + itos( *varIter );
	varIter++;
	while( varIter != edges->end() )
	{
		problemStr += "\t+\tX" + itos( *varIter );
		m_constraints += "\t-\tX" + itos( *varIter );
		varIter++;
	}
	problemStr += "\n";
	m_constraints += "\t=\t0\n";
}

void ContigGraphAnalyzer::GenAimFunction( Gap *gap, string &problemStr )
{
	set<Contig*> *graph = gap->GetSubgraph();
	problemStr += "\\Aim:\n\\";
	for( set<Contig*>::iterator conIter = graph->begin(); conIter != graph->end(); conIter++ )
	{
		list<int> *inEdgeVars = (*conIter)->GetNewInEdgeVar();
		for( list<int>::iterator varIter = inEdgeVars->begin(); varIter != inEdgeVars->end(); varIter++ )
		{
			problemStr += "X" + itos( *varIter ) + "*" + itos( (*conIter)->GetLength() - Configure::OVERLAP ) + "\t+\t";

			// generate quadratic coefficients
			if( m_aimStringLinear == "" )
				m_aimStringLinear = "-";
			else
				m_aimStringLinear += " - ";
			m_aimStringLinear += itos( (gap->GetGapSize() + Configure::OVERLAP) * 2 * ((*conIter)->GetLength() - Configure::OVERLAP ) ) + " X" + itos( *varIter );
			AddAimLinearCoef( -(gap->GetGapSize() + Configure::OVERLAP) * 2 * ((*conIter)->GetLength() - Configure::OVERLAP ), *varIter );
			int len = (*conIter)->GetLength() - Configure::OVERLAP;
			m_aimStringQuad += itos( len * len * 2 ) + " X" + itos( *varIter ) + " ^ 2 ";
			AddAimQuadCoef( len, *varIter, *varIter );
			for( set<Contig*>::iterator nextConIter = graph->begin(); nextConIter != graph->end(); nextConIter++ )
			{
				list<int> *inEdgeVarsNext = (*nextConIter)->GetNewInEdgeVar();
				for( list<int>::iterator varIterNext = inEdgeVarsNext->begin(); varIterNext != inEdgeVarsNext->end(); varIterNext++ )
				{
					if( *varIterNext != *varIter )
					{ // different variable
						if( *varIter < *varIterNext ){
							if( m_aimStringQuad != "" )
								m_aimStringQuad += " + ";
							m_aimStringQuad += itos( ((*conIter)->GetLength() - Configure::OVERLAP) * ((*nextConIter)->GetLength() - Configure::OVERLAP) * 4 ) + " X" + itos( *varIter ) 
								+ " * X" + itos( *varIterNext );
							AddAimQuadCoef( ((*conIter)->GetLength() - Configure::OVERLAP) * ((*nextConIter)->GetLength() - Configure::OVERLAP) * 4, *varIter, *varIterNext );
						}
					}
				}
			}
		}
	}

	problemStr += "\t-\t" + itos( gap->GetGapSize() + Configure::OVERLAP ) + "\n";
}

void ContigGraphAnalyzer::AddAimQuadCoef( int value, int x, int y )
{
	// always keep x <= y
	if( x > y )
	{ // swap x, y
		int temp = x;
		x = y;
		y = temp;
	}

	while( x + 1 > m_aimQuadValue->size() )
	{
		m_aimQuadValue->push_back( new vector<int> );
	}

	vector<int> *xVector = m_aimQuadValue->at( x );
	while( y + 1 > xVector->size() )
	{
		xVector->push_back( 0 );
	}

	(*xVector)[ y ] += value;
}

void ContigGraphAnalyzer::AddAimLinearCoef( int value, int id )
{
	while( id + 1 > m_aimLinearValue->size() )
	{
		m_aimLinearValue->push_back( 0 );
	}

	(*m_aimLinearValue)[ id ] += value;
}

void ContigGraphAnalyzer::CreateVariable( Contig *c, map<int, pair<Contig*, int> > *edges, int type, set<Contig*> *subgraph, vector<int> *variables, string &problemStr, vector<int> *scafVar )
{
	int tailType, headType;
	if( type == LEFT )
		tailType = IN;
	else
		tailType = OUT;

	bool inScaf = false;
	if( c->IfInScaffold() )
		inScaf = true;
	
	// check each edge
	for( map<int, pair<Contig*, int> >::iterator iter = edges->begin(); iter != edges->end(); iter++ ){
		if( c->GetID() > iter->second.first->GetID() )
			continue;

		if( type == RIGHT && c->GetID() == abs( iter->first ) )
			continue;

		Contig *otherContig = iter->second.first;

		// do not output the contigs which is not in this subgraph
		if( subgraph->find( otherContig ) == subgraph->end() )
			continue;

		// do not output contigs with low coverage
		if( otherContig->GetCov() < m_covThreshold )
			continue;

		// if it is within the subgraph
			// check the edge type for other contig
		if( type == RIGHT ){
			if( iter->first > 0 )
				headType = IN;
			else
				headType = OUT;
		}
		else{
			if( iter->first > 0 )
				headType = OUT;
			else
				headType = IN;
		}

		// create a variable and save to both contigs
		if( tailType == IN )
			c->AddNewInEdgeVar( variables->size() );
		else
			c->AddNewOutEdgeVar( variables->size() );

		if( headType == IN )
			otherContig->AddNewInEdgeVar( variables->size() );
		else
			otherContig->AddNewOutEdgeVar( variables->size() );


		problemStr += "\\Description:\t";
		problemStr += "X" + itos( variables->size() ) + "\t";
		if( tailType == IN )
			problemStr += "-";
		problemStr += c->GetName() + "\t";
		if( headType == OUT )
			problemStr += "-";
		problemStr += otherContig->GetName() + "\n";

		m_variableString += " X" + itos( variables->size() );
		
		// check if related with start or end contig
		if( inScaf )
			scafVar->push_back( variables->size() );

		variables->push_back( variables->size() );
	}
}

// add the punishment function
// c*[sum(x)(x1+...+xn) - 2(x1+...+xn)]
void ContigGraphAnalyzer::AddPunishFunction( vector<int> *scafVar, vector<int> *var )
{
	int c = 99;
	for( int i = 1; i < var->size(); i++ )
	{
		AddAimLinearCoef( 2*c, i );
	}

	for( int i = 0; i < scafVar->size(); i++ )
	{
		for( int j = 1; j < var->size(); j++ )
		{
			this->AddAimQuadCoef( -2*c, scafVar->at( i ), j );
		}
	}
}

// the pipeline of analyzing contig graph
int main(int argc, char *argv[] )
{
	ContigGraphAnalyzer *m_analyzer = new ContigGraphAnalyzer();
	//m_analyzer->m_fileType = ABYSS;					// the type of input file, soap or velvet
	m_analyzer->m_fileType = VELVET;
	//m_analyzer->m_fileType = SOAP;

	/*if( argc != 3 && argc != 4 )
	{
		cout<<"Usage:\n";
		cout<<"Please specify the SOAPdenovo related files as follows:\n";
		cout<<"If the threshold of coverage is not specified, all the contigs will be considered\n";
		cout<<"\tDotGenerater <input data prefix> <output folder> (<coverage threshold>)\n";
		return 0;
	}

	string m_inputPrefix = argv[ 1 ];

	m_analyzer->m_outputPath = argv[ 2 ];
	mkdir (m_analyzer->m_outputPath.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
	if( m_analyzer->m_outputPath.at( m_analyzer->m_outputPath.length() - 1 ) != '/' )
		m_analyzer->m_outputPath += "/";

	string m_graphFileName, m_contigFileName;
	if( m_analyzer->m_fileType == SOAP ){
		m_graphFileName = m_inputPrefix + ".updated.edge";
		m_contigFileName = m_inputPrefix + ".contig";
		// read kmer size
		if( m_analyzer->GetKmerOfSOAP( m_inputPrefix ) == -1 )
			return -1; 
	}
	
	if( argc == 4 )
		m_analyzer->SetCovThreshold( atof( argv[ 3 ] ) );
	else
		m_analyzer->SetCovThreshold( 0 );*/

	string m_graphFileName = "ecoli-contigs.dot";
	string m_contigFileName = "ecoli-contigs.fa";
	m_analyzer->SetCovThreshold( 0 );
	m_analyzer->m_outputPath = "./";

	// step 1: Read config file
	string configFileName;
	if( argc == 2 ){
		configFileName = string( argv[ 1 ] );
		ConfigReader myConfigReader;
		cout<<"Step 1: read configuration file"<<endl;
		if( myConfigReader.ReadConfigFile( configFileName ) == -1 )
		{
			cout<<"Reading Configuration file unsuccessfully!\n";
			return -1;
		}
	}

	// step 1: read graph file to get the information of contigs
	cout<<"Step 2: read graph file to get the information of contigs and edges"<<endl;
	if( Configure::GRAPH_TYPE == VELVET ){
		if( m_analyzer->ReadVelvetGraphFile( Configure::DATA_DIRECTORY + "LastGraph" ) == -1 ){
			cout<<"\tError reading graph file"<<endl;
			return -1;
		}

		// step 1.5: read contig file to get the contig name
		cout<<"\treading contig fasta file...\n";
		if( m_analyzer->ReadContigFasta( Configure::DATA_DIRECTORY + "contigs.fa" ) == -1 ){
			cout<<"\tError reading contig fasta file\n";
			return -1;
		}
	}
	else if( Configure::GRAPH_TYPE == SOAP ){
		if( m_analyzer->ReadSoapFile( m_graphFileName, m_contigFileName ) == -1 ) {
				//"/mnt/ScratchPool/bertrandd/WORM/ASSEMBLY/WOA001/RAW/WORM.updated.edge", "WORM.contig" ) == -1 ){  
			//     "/home/bertrandd/citrus/quake_filtring/soapdenovo/k_41/CITRUS.updated.edge", "CITRUS.contig" "/mnt/ScratchPool/bertrandd/WORM/ASSEMBLY/WOA001/RAW/WORM.updated.edge", "WORM.contig" "test.edge", "test.contig"
			cout<<"\tError analyzing SOAPdenovo file\n";
			return -1;
		}
	}
	else if( m_analyzer->m_fileType == ABYSS ){
		if( m_analyzer->ReadAbyssFile( m_graphFileName, m_contigFileName ) == -1 ) {
			cout<<"\tError analyzing ABySS file\n";
			return -1;
		}
	}

	// calcuate average coverage and set the copy number for each contig
	m_analyzer->CalculateAvgCov();

	// step 2: output dot file
	// cout<<"output graph file\n";
	//cout<<"Step 2: output dot file\n";
	//m_analyzer->PrintDotFile();

	// print graph around vector
	//m_analyzer->PrintGraphAroundVector();
	
	// step 5: analyze the subgraph of each gap
	cout<<"Step 5: Analyze subgraph of each gap and find the connectivity\n";
	m_analyzer->AnalyzeGaps( Configure::SCAFFOLDS_FILE );

	delete m_analyzer;
}
